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Abstract 

This reports deals with three major topics related to TSS dynamics: 

a) The update of the SAO rotational dynamics computer code. The pro- 
gram is now suitable to deal with inclined orbitB. The output has been 
also modified in order to show the satellite Euler angles referred to the 
rotating orbital frame. 

b) The development of the three-dimensional high resolution computer 
program SLACK3 . The code simulates the three-dimensional dynamics of a 
tether going slack taking into account the effect produced by boom rota- 
tions. Preliminary simulations on the three-dimensional dynamics of a re- 
coiling slack tether are shown in this report. 

c) The development of a program to evaluate the electric potential 
around a severed tether immersed in a plastr.a . The potential is computed on 
a three-dimensional grid axially symmetric with respect to the tether lon- 
gitudinal axis. The electric potential variations due to the plasma are 
presently under Investigation. 
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1.0 INTRODUCTION 

This Is the second quarterly report submitted by SAO under contract 
NASS-36160, "The Investigation of Tethered Satellite System Dynamics," Dr. 
Enrico Lorenzini, PI, and covers the period from 15 November 1984 through 
16 February 1985. 


2.0 TECHNICAL ACTIVITY DURING REPORTING PERIOD AND PROGRAM STATUS 

2.1 Tss Rotational Dynamics Analysis 

2.1.1 Introductory Remarks - 

The SKYHOOK computer code was updated with a preliminary rotational 
dynamics package in 1982. This study is documented in the report "Study of 
Tethered Satellite Active Attitude Control," October 1982, NASA Contract 
NAS8-33691. The satellite was modeled as a rigid body and the tension 
transmitted by the tether was the only external perturbation acting upon 
the satellite. The three Euler rigid body equations were added to SKYHOOK 
along with the nine equations relating the body-axis angular velocities to 
the derivatives of the direction cosines. Once the direction cosines are 
obtained, the Euler angles are computed. This method has the major draw- 
back of integrating nine additional equations that can be reduced to three. 

The body-axis angular /elocitles , Wi, Wj, W 3 of the Euler rigid body 
equations can be converted to Euler angular velocities and then Integrated 
a second time to get the Euler angles. The direct integration of the Euler 
angular velocities will therefore increase the computer code efficiency and 
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reduce the CPU time. The only drawback related to this approach is that 
the equations become singular when the roll angle approach 0°, that is to 
say (with the SKYHOOK reference frame conventions) when the vertical body 
axis of the satellite is perpendicular to the equatorial plane. This con- 
dition however is a limit condition that can be encountered only in special 
situations. Both integration options will be kept available in the SKYHOOK 
computer code to have the greatest flexibility. 

The other major modification is tne implementation of a new force 
package and of a new torque package that models the effect of thruster 
activation. In-plane, out-of-plane and in-line thrusters provide, in an 
ideal situation, forces applied to the satellite center-o f-mass (c.m.) . In 
the actual situation they also produce spurious torques because of 
misalignments. These torques must appear on the right hand side of the 
rigid body Euler equations. However, this doesn't mean that a thruster 
activation produces a satellite rotation only through misalignment. Even 
If a lateral force acting upon the satellite is directed exactly toward the 
satellite c.m., rotations are still produced by the different moments of 
inertia of the satellite and the tether. This effect is modeled in SKYHOOK 
so that the tether bends differently depending upon the point of applica- 
tion of the external force and the mass distribution of the system. 

The update of the output is also necessary in SKYHOOK. Presently the 
computer code plots the satellite attitude with respect to the inertia] 
reference frame. Plots cannot be interpreted immediately especially if the 
orbit is not equatorial. Reference frame transformations are necessary to 
convert the inertial Euler angles to the rotating reference frame Euler 
angles. In this way the satellite attitude is referred to the local verti- 
cal and the local horizon and the output has an immediate interpretation. 
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This way of dealing with inertial Euler angle has made it difficult to set 
up the initial conditions for SKYHOOK. For this reason, the rotational 
dynamics SKYHOOK has been used, in the past, for equatorial orbit only (in 
this case the problem is greatly simplified). The matrix transformation 
implemented to modify the output can also be used for setting up initial 
conditions in order to generalize the program to any type of orbit. All 
these modifications are in progress. Some of them will be fully operative 
in our computer code in the near future. 


2.1.2 Matrix Transformation From Inertial Reference Frame To Orbiting Ref- 
erence Frame - 

The attitude dynamic equations integrated by SKYHOOK make use of the 
direction cosine matrix. From that matrix we derive the Euler angles with 
respect to the inertial reference frame (x,y,z) (see Figure 2.1.1). If we 
want to refer the satellite attitude to the orbiting reference frame (x* , 
yl , 7.1) we must derive the rotation matrix that relates the orbiting refer- 
ence frame to the inertial reference frame. If we call [Roi] the rotation 
matrix that transform orbiting axis into inertial axis and | Rbi | the rota- 
tion matrix that transform satellite body axis into inertial axis we have: 


Xi 


• 

X B 

Yi 

= [Rbi 3 

ya 

. 21 . 


_Z B . 


( 2 . 1 . 1 ) 
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and 


Xi 


M 

yi 

11 

r — 1 

50 

0 

M 

V* 

Zl 


Z> 


so that 


'xf' 


' 

Xo 

Y' 

= [Roi ] 1 [Rax] 

yu 

zf 




where 


( 2 . 1 . 2 ) 


(2.1.3) 


[Roi] T - [Rio] (2.1.4) 

is the transposed of the matrix [Roj] while (xb, ye, z B ) is the satellite 
body-axis reference frame as shown in Figure 2.1.1. The matrix [Roi] is 
derived through three sequential rotations of the inertial reference frame 
(xi, yj , Zj) in order to transform it into the orbiting reference frame 
(xf , y> , 2.1). The first rotation is around the axis z r and defines the 
line of nodes of the orbit so that the angle is the argument of nodes, 

the second one is around the transformed axis y z f and defines the orbit 
inclination i and the third one is around the transformed axis z z l and 
defines the argument of latitude 7 . Because of the orientation of the 
orbiting reference frame adopted in SKYHOOK and depicted in Figure 2.1.1 
the inertial reference frame needs to be further rotated 180° around the 
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axis x> . This extra rotation implies a change in sign of the second and 
third row of the rotation matrix so that: 



"l 

0 

0' 


cosy 

siny 
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and finally: 


[Rio] 


cosycos/9- 

cosysin/3+ 

sinysini 

sinysin^cosi 

sinycos^cosi 


sinycos/J* 

sin^siny- 

-cosysini 

cosycosiain^? 

cosycos/JcoBi 


-sinisin/J 

sinicos/? 

-cosi 


( 2 . 1 . 6 ) 


The rotation matrix between the satellite body reference frame and the 
inertial reference frame was derived in the Interim Report "Study of Teth- 
ered Satellite Active Attitude Control," G. Colombo, October 1982, Contract 
NAS8-33691, The reverse transformation matrix, from body-axis to inertial 
axis is as follows: 


. '■ir - ■ 
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cos^cos <(>- -sin^cos^- sintfsin^ 

cos^Bln^sinV* cos0Bin$cos^ 


[RbiJ - |co8^oin^+ 

j costfcos^sin^ 

sin^sintf 


-sin^sin^t -sinflcos^ 

COS0COS^COS\& 


cos^sin0 


cos 0 


(2.1.7) 


In equation (2.1.7) ip is the pitch angle, 0 is the roll and <P is the yaw 
angle of the satellite. All of them are referred to the inertial reference 
frame so that they give the inertial attitude of the satellite. 

The matrix to transform the body reference frame into the rotating 
reference frame is therefore given by: 


[Rdo] = [Rio] [Rdi] 


( 2 . 1 . 8 ) 


where the matrix [Rio] and [Ro:] are given by expression (2.1.6) and 
(2.1.7) respectively. The matrix multiplication is not performed analyti- 
cally since it in more convenient to implement it on the computer. Once 
the matrix [R D o] is obtained ws can derive the Euler angles with respect to 
the orbiting reference frame as follows: 


[Rbo 


ai 

a J 

Sj 

Pi 

Pi 

Pi 

'll 

72 

73 


and 


(2.1.9) 
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tan <f>l = - — 

fi s 

tan V" 1 = — 1 

72 

tan Oi =* \/'i) a +'u*/ 7s (2.1.10) 

The last relation has an unavoidable! sign ambiguity. The sign is chosen 
positive for convenience. 

The matrix [Rbo] can also be written starting from the Euler angles 
relative to the orbiting reference frame (tpl , 6> t <f>>) . Referring to equa- 
tion (2.1.9) the elements of the matrix can be written as follows: 

a i = coBtpi coB<f>> - cob9 I Blmj/t 
a 3 = -sinV ,/ cos^ - cosO I sin ^* coBipl 
a 3 ~ sin0 1 sin#* 

/3i = cos^sin^f + cobO 1 coBif>i si.n^/1 
f} 3 = -sin^'sin^f + cos0 *cos^* cos^' 

— -sinfMcos^f 
7i = Bin^'sinfl* 

7a = cos^fsin## 

7 a = cosfl 1 (2.1. 11) 


Because of equation (2.1.8) the matrix of the body- axis direction cosines 
in inertial axis is given by: 


[Rbi] = [Rio] t [Rbo] 


( 2 . 1 . 12 ) 












In equation (2.1.13) (wjf , oi 3 * , w 3 *) are the body-axis angular velocity 
components as measured from an observer moving with the rotating reference 
frame; («i , Wj , w 3 ) are the inertial angular velocity components in body- 
axis; (11*, Dy, 0,) are the inertial angular rates of the rotating refer- 
ence frame in inertial axis while the matrix [Rbi3 t = [Rib] transforms the 
components from inertial axis to body-axis. The vector (u>i, u 3 , a>s) has a 
special significance because its components are the quantities measured by 
a gyro package strapped down to the satellite. (u>i/ , , w 3 i) are similar 

components but they are obtained from the inertial vector (uj_, w 3 , w 3 ) 

after filtering out the orbital motion of the system. 

All these matrix transformations allows the simulation program to 
switch from the inertial reference frame representation to the orbiting 

I: 

.. _ „-*r ••c,-?-’-:- — ' 

.ih — - — ■ — — -- 






•Z'2^ 


Page 12 

reference frame representation and vice versa. The outcome is to present 
the output in a physically meaningful form and to have the possibility of 
setting up initial conditions in SKYHOOK for a generic orbit. 



Figure 2.1.1 Satellite Rotational Dynamics. 

System Reference Frames. 
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2.1.3 Relations Between The Euler Angular Velocities And The Body-Axis 
Angular Velocities - 

In the present version of SKYHOOK the body-axis angular velocities are 
derived from the integration of Euler equations of the satellite. The 
body-axis angular velocities are then related to the derivatives of the 
direction cosines . After a second integration the direction cosine matrix 
is obtained and the inertial Euler angles are derived by means of equations 
like (2.1.10). This process can be strongly simplified as follows: a) 
after the first integration of the Euler equations, the body-axis angular 
velocities can be analytically transformed to Euler angular velocities, b) 
inertial Euler angles can then be directly obtained with a second integra- 
tion . 


The body-axis angular velocities are related to Euler angular velocl 
ties as follows: 


Wl 


sin#sint# coBrf/ 0 



Wj 

= 

sin^cos^ -sin^ 0 



_«3_ 


cos# 0 1 




(2.1.14) 


If the transformation matrix is reversed (note that this is not an orthogo- 
nal matrix) we get: 
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cos\fr/sin 6 0 
O 

-coBipcotgO 1 

The problem associated with equation (2.1.15) is that it becomes singular 
when the inertial roll angle 0 = 0° that is to say when the satellite 
longitudinal axis is perpendicular to the equatorial plane. This means, 
for example, that in the case of an equatorial orbit, equation (2.1.15) 
becomes singular when the satalllte is at 90 degree with respect to the 
tether (because of the reference frame conventions in SKYHOOK 6 = 0° in 
this case) . As the orbital inclination Increases the margin on the iner- 
tial roll angle that makes equation (2.1.15) singular is reduced. Far from 
the polar orbit, however, there are no problems in running SKYHOOK with 
this more efficient procedure. 

The line that we are following is to have both options available in 
the computer code in order to use the novel and more efficient approach 
when possible. The standard option will be used for the special cases when 
the satellite attitude is near a singularity. 


Wl 

w 2 

«3 


(2.1.15) 
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sin^/sin0 
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-sin^cotgtf 
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2.1.4 Torque And Force Package Update - 

The external torques acting upon the satellite are as follows: 

a. The yaw torque provided by the yaw control thrusters 

b. The spurious torques due to misalignments of the out-of-plane, 
in-plane and in-line thrusters 

c. The torques due to the passive dampers of the satellite 

d. The torque due to the tether tension 

Since, in our computer code, the satellite attitude dynamics is modeled, by 
Euler equations, all the external torques must be expressed in satellite 
body-axis. The Euler equations of the satellite are given by: 

IiWx - W3W3 (I 3 -I 3 ) 

IjW 2 - U)U) 1 (Ij-Ii) 

Ijw 3 - WjWj (Ix-Ij) 

In (2.1.15), the vector N T{m is the torque due to the tether tension (it is 
already modeled in the computer code) , the vector Nyi Is the torque due to 
the yaw thrusters, the vector N„i S is the torque due to thruster misalign- 
ments and N Bm is the one due to the dampers on board the satellite. 

Iii Iji I 3 ate the satellite moments of inertia and Wi, u>j, w 3 are the 
components of the satellite rotational velocity in body-axis. The torque 
Nyt has one component only, along the vertical axis of the satellite, so 


= N Ten + N n + N mib + N d . 0 (2.1.16) 


that: 
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0 



(2.1.17) 


Njyr must be modeled according to the latest version of the AERITALIA con- 
trol logic. In general it is a pulse modulated signal, presently with an 
amplitude of 1 Nm, Pulses are modulated according to the yaw angle <f> and 
the yaw angular rate <j> in order to stabilize the yaw attitude of the satel- 
lite or to maintain a constant spin rate. 


The torque N M ib can be itemized in terms of in-plane, out-of-plane and 
in-line thrusters as follows: 
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( 2 . 1 . 18 ) 


In equation (2.1.18) F IP , E 0 p and Fi L are the thrust level of the in -plane, 
out-of-plane and in-line thruster respectively. The vectors bx , b a and b 3 
are the arms of the thrust vectors with respect to the satellite c.m. 
Since the satellite is spherical and the thrusters are located on the outer 
shell, the arm vector can be expressed in terms of angular misalignments as 


bk = 


follows : 


R 


sin£*k,x 

sina k , 2 

sina k ,3 


(2.1.19) 
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where : 

k = (1,2,3) and a kjj = 0 if k = ] (2.1.20) 

a ktJ ?! 0 if k ?! j 

a kt j are the cartesian components of the angular misalignment between the 
actual thruster longitudinal axis and its theoretical position. In the 
case of an in-plane thruster, for example, the overall angular misalignment 
is given by: 


| oti | = sin " 1 [sin 2 (a 1>a ) + sin 2 ( 01 , 3 ) ] ( 2 . 1 . 21 ) 

where x B = 1 , yn = 2 and z B = 3 ; x B , y B and z B , being the satellite body 
axis . 


The thrust levels Fjp, F 0 p and Fn, must be modeled according to the 
selected control logic. In particular the in-line thruster has a two- 
staged thrust level and it is switched an, in a continuous operation mode, 
every time that the tension fetlls below 2N as follows: 

F il = IN if IN < T < 2N 
Fil = 2N if T < IN (2.1.22) 

where T is the tether tension . 

The most efficient way to activate the out-of-plane thrusters is when 
the satellite crosses the local vertical so that its velocity is maximum. 
To further improve the efficiency it is convenient to let the out-of-plane 
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oscillation build up to 20° or 30° and then fire the thruster at the next 
crossing. The firing time is given by: 

Atop = {®bh + ^ n't) Vopoax/Fop ~ {m»a + mt) (2fl B bbx £)/Fop (2.1.23) 

In equation (2.1.23) m (I is the satellite mass, m t is the tether mass (that 
is negligible for short tether length) , V 0Pmx is the out-of-plane velocity 
of the satellite at the crossing of the local vertical, B nxx is the ampli- 
tude of the out-of-plane angle of the tether before thruster activation. 
All the other terms are self explanatory. It is interesting to notice that 
the satellite velocity V 0Pbbx is fairly constant with the tether length, 
In fact from the out-of-plane rotational energy conservation (it holds true 
in first approximation because the damping is very low) we get: 

Vopoax = 2nB o £ f ,/[(m 0fl + | nC)/( m aa + | fit 0 ) ] V* (2.1.24) 

In equation (2.1.24) the index "o" stands for values computed at the begin- 
ning of retrieval and (i is the linear density of the tether. For tethers 
few kilometers long, equation (2,1.24) is essentially independent from the 
tether length l . 

The in-plane thruster control logic is theoretically similar to the 
out-of-plane one. However, the in-plane control is more complicate because 
the tether oscillation has a variable bias dependent upon the Coriolis 
force and the air drag. The tiring time will have an expression similar to 
equation (2.1,23), as follows: 

Atjp = (m aa + —mt) Vi Pmax /Fip — (m 00 + ^ O/Bpp (2.1.25) 




< 3 * 
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This time, however, all the angular quantities must be computed from the 
in-plane bias and the in-plane thrusters must be activated when the satel- 
lite crosses it. The evaluation of the in-plane oscillation bias could be 
performed by software processing of the radar observations (if this is at 
all feasible). A simplified analytical expression for it, valid for an 
exponential retrieval control law, is as follows: 



sin 1 


[■ 


2 

t c ficos (At) 


1 s 

2 m, B 


CdPa 


(R.+H) » 

t J 


(2.1.26) 


Equation (2.1.26) can be further simplified, for small angle, by assuming 
cos (Ab) a 1 so that it becomes very easy to compute the angular bias. In 
equation (2.1.26) t c is the retrieval time constant, C D is the drag coeffi- 
cient of the satellite, S its cross section (the tether cross section is 
neglected because the bias angle builds up whan the tether is short) , p A is 
the atmospheric density, R e the earth radius and H the orbital altitude. 

The software Implementation of the torque and force package in our 
computer code is presently underway. The torque N Ten due to the tether 
tension is already available. The torque N Danp due to the satellite atti- 
tude dampers has been examined only in a simplified way. The implementa- 
tion and design of roll and pitch passive dampers has been recently studied 
by the Italian prime contractor. We will contact AERITALIA in the near 
future to understand what the present status is in order to model the 
damping torques accordingly. 
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2.1.5 Current Status Of Rotational Dynamics Code - 

Under NASA Contract NAS8-32199, a version of the DUMBEL program was 
created which includes the rotational dynamics of the subsatellite. This 
code has a model for a damping torque proportional to the rate of change of 
the angle between the wire and the vector from the center of mass of the 
subsatellite to the attachment point of the wire. The program is docu- 
mented in Appendix III of the Final Report "Study of the Dynamics of a 
Tethered Satellite System (SKYHOOK)," Kalaghan et al., March 1978. 

Under NASA Contract NAS8- 33691, the code developed in DUMBEL, but 
without the damping model, was implemented in the SKYHOOK computer program. 
This work is described in the Interim Report, "Study of Tethered Satellite 
Active Attitude Control," G. Colombo, October 1982. The code was tested 
with a rudimentary thruster model that applied a constant torque about any 
of the three body axes. A special test case was also devised in which the 
subsatellite precesses under the torque of the wire at the same rate as the 
orbital angular velocity. In this case, the subsatellite maintains a fixed 
orientation between the symmetry axis and the rotating orbital reference 
frame. In order to show the orientation in the orbital frame, an equato- 
rial orbit was run and the direction cosine matrix was rotated about the z- 
axis by the orbital angle in order to derive Euler angles in the rotating 
reference frame. The code was not developed to the point of being able to 
give Euler angles with respect to an arbitrary orbital reference frame. 
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2.1,6 Implementation Of The Transformation To The Rotating Orbital Refer- 
ence Frame - 

In order to perform simulations for arbitrary orbital elements and 
obtain meaningful displays of the subsatellite rotation it is necessary to 
be able to set up initial conditions in a rotating orbital reference frame 
and be able to display the output orientation angles in that coordinate 
system. In the method of direction cosines, the coordinates of a point 
with respect to the body axes are given as a function of the inertial 
coordinates x by the equation 


x' -Ax 


(2.1.27) 


where 


<*1 

P i 

'll 

aj 

Pi 

li 

a 3 

Pi 

li 


(2.1.28) 


The matrix A can be expressed as the product of three matrices representing 
rotations by the Euler angles. That is, 

A = B(^)C(0)D(^) (2.1.29) 


where 


BW = 


coeip sLnip 0 
-sin ip cosip 0 


0 


0 


1 



In the case of an equatorial orbit with the Shuttle initially on the 
x-axis, the orbital reference frame initially coincides with the Inertial 
reference frame so that equation (2.1.29) gives the initial conditions. In 
general the transformation from the inertial frame to the rotating orbital 
reference frame can be given as: 


l 

y 

U 


i 


E = B( 7 )C(i)D(/?) (2.1.30) 

where j3 is the argument of the ascending mode, i is the orbital inclina- 
tion, and 7 is the argument of latitude (the angle from the node to the 
Shuttle in the plane of the orbit). Equation (2.1.30) has the same form as 
equation (2.1.29). The transformation E gives the coordinates of a point 
in the orbital reference frame as a function of the coordinates in the 
inertial frame. For arbitrary orbital elements the transformation from 


F = GE 




(2.1.31) 
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and 


G = B(^)C(0')D(tf') (2.1.32) 

The angles 6> , <f>< , and if>* give the relationship between the orbital frame 
and the body axis frame. Equation (2.1.31) giving E as a function of six 
rotations can be used to set up the initial direction cosine matrix for 
arbitrary orbital elements and arbitrary orientation of the body axes with 
respect to the orbital reference frame. 

The tension in the wire is computed using the distance from the Shut- 
tle to the attachment point of the wire. In order to set up equilibrium 

initial conditions, the positions of the attachment point and the center of 
mass of the subsatellite must be known in inertial coordinates. If the 
coordinates of the attachment point in the body axis frame are given by the 
vector r> , the coordinates in inertial space are given by 

r = F T r' (2.1.33) 

where E T is the transpose of the matrix F. In the DUMBEL program the 
vector r is subtracted from the position pf the attachment point to compute 
the initial position of the center of mass of the subsatellite. 

The rotational motion of the subsatellite is obtained by integrating 
vs time the nine elements of the matrix F and the three components of the 
angular velocity along the body axes. The angular quantities are not inte- 
grated directly. At any instant of time the nine direction cosines can be 
Interpreted either as the results of three rotations as in the case of the 
matrix A, or as six rotations as in the case of matrix F. A convenient 


*2r •*£.&< 
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form of output would be to give the three angles, 0>, , and tj>l which 
define the transformation G from the orbital frame to tody axes. Since F = 
GE we can obtain the matrix G by post multiplying by E T to give 

EE 1 = GEE t = G (2.1.34) 

We can then derive the angles Oi, <f>i , and ipi from the matrix G. A method 
of obtaining these angles is derived in Section 2.1.2 of this report. The 
results of the derivation are 


tan^f = -a 3 /f} 3 

tan it>f = 7i/7a 

tan0* = \/7i a + 72 s /7s 


(2 .1 .35) 


where a and 7 are direction cosines as defined for the matrix A. The 
explicit representation of the direction cosines in terms of Euler angles 
is given in equations (2.1.11) of the above cited Section. 


2.1.7 Rotational Dynamics Co^ j First Level Update - 

As a first step in the present study, a few of the cases used to test 
the program have been rerun and the output checked against the original 
test runs. The original version of the program uses explicit expressions 
for the elements of the matrix in equation (2.1.28) . The matrix F involves 
six rotations, which would be very cumbersome to do analytically. There- 
fore, a set of general purpose subroutines have been written to manipulate 
3x3 matrices numerically. Subroutine ROTMAT (IAXIS , ANGLE ,A) sets up a 
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matrix A to perform a rotation about axis number "I AXIS" through an angle 
"ANGLE." The explicit form of the matrices for each axis with ANGLE = ® are 


1 = 


10 0 
0 cos® sin® 

0 -sin® cos® 


2 = 


cos® 

0 

sin® 


0 

1 

0 


-sin® 

0 

cos® 


3 = 


cos® 

-sin® 


sin® 0 
cos® 0 
0 1 


Subroutine MATMUL (A,B|C) performs the matrix multiplication C — A x 
B. Subroutine STATVEC (THETA, PUI , PSI , A) creates the direction cosine 
matrix A = B (l®) C (®) D ( 4 >) . This has been tested against the original subrou- 
tine which uses analytic expressions for the nine direction cosines and 
gives the same results. Subroutine TRNSPOS (A, AT) performs the operation 
AT = A T . Subroutine VECMUL(A. V, VP) performs the operation VP = A x V 
where V and VP are vectors and A is a 3x3 matrix. The original version of 
subroutine ROTSTAT in program DUMBEL has been replaced by a more general 
version ROTSTAT (AN,RI , AL , THETA , PHI ,PSI ,DTHETA,DPHI ,DPSI ,Y,A,DIRC0S) , where 
AN is the argument of the node, RI is the orbital inclination, AL Is the 
argument of latitude, THETA, PHI and PSI are the Euler angles with respect 
to the rotating orbital coordinate frame. DTHETA, DPHI , and DPSI are the 
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rates of change of the angles in inertial space but referred to the axes of 
the rotating orbital reference frame. The angular velocity in inertial 
space is the angular velocity in the rotating frame plus the orbital angu- 
lar velocity. The way the Euler angles are defined, is parallel to the 
orbital angular velocity, so that the angular velocity in inertial space 
can be obtained by adding the orbital angular velocity (1 to the <pl which is 
the rate of change of <f>> with respect to the rotating orbital frame. The 
components of the angular velocity along the body axes can be obtained 
using the equations 


Wj. = ^sin0 fsint/h + $cosifii 

(1)2 — <J>sLnO i cosip 1 - $sin^< 

a) 3 = ^cos<M + rj) 

where 6, /}>, and (j) are the rates in inertial space. 

Subroutine ROTSTAT creates the matrices 

E = B (AL) C (RI ) D (AN) 


and 


G = B(f )C(^)D(^') 

using subroutine STATVEC and then performs the operation F = GxE . The 
matrix E is returned as DIACOS , and placed in A, the rotational part of the 
state vector. (The positional part, Y, of the state vector is not used in 
this version) . The position of the attachment point of the wire is com- 
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Y'» 


tv 


puted by the equation 


R t - F^b 


and returned via common . 

In order to display the rotation of the subsatellite in the rotating 
orbital reference frame, we must have the instantaneous orbital parameters. 
These can be obtained irom the position P and velocity V. The unit vector 
normal to the orbit n is in the direction PxV. The orbital inclination is 
cos' 1 n s . The argument of the node is tan" 1 n y /n*. If we call x a unit 
vector along the node and y a unit vector in the direction n x x, we can 
get the sine and cosine of the argument of latitude from the components of 
P along y and x. Subroutine ELEM (P , V,AN ,RI ,AL) computes the argument of 
the node AN, the inclination RI , and argument of latitude AL from the 
position P and velocity V. 

A routine called ANGROT has been written previously to rotate the 
direction cosine matrix about the z-axls and compute Euler angles from the 
rotated matrix. This routine has been used for equatorial orbits to give 
the orientation of the subsatellite in a rotating orbital reference frame. 

A new version ANGROT (AN ,RI ,AL ,A,ANG) has been written where AN is the 
argument of the node, RI is the inclination, and AL is the argument of 
latitude, A is the instantaneous direction cosine matrix called E in the 
general case, and ANG is a three element array containing the Euler angles 
with respect to the orbital reference frame. Subroutine ANGROT performs 
the operation 


G = EE T 
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I 

using subroutines STATVEC, TRNSPOR, and MATMUL. The vector ANG is then 
| obtained from G. The calls to ELEM and ANGROT are coordinated by subrou- 

tine ROTCDORD (Z,A,ANG) , where Z is the positional part of the state vec- 
tor, A is the matrix of direction cosines and ANG is the Euler angles in 

1 

} the rotating orbital reference frame. 



2.1.8 Rotational Dynamics Computer Simulations - 

Various test cases have been run with the new software described in 
Section 2.1.6. In the first case the system is in an equatorial orbit with 
no initial rotation or angular velocity in inertial space. Angular damping 
is included in the model. In the original version of the program, the 
rotation of the wire at the orbital angular velocity causes the subsatel- 
lite to acquire the orbital angular velocity. Measured in the orbital 
reference frame, the subsatellite initially rotates away from equilibrium 
and then returns exponentially to equilibrium. The new software produces 
the same results and agrees numerically with the previous results. 

The second test case is the same as the first but with an orbital 
inclination of 28°. This run appears to have integrated properly but the 
output angles were disrupted by a singularity in the angular representa- 
tion. When the angle 6 1 is zero, one cannot distinguish between 4>> and *j)t . 
In this case the program arbitrary sets <f>> =0 and attributes all the 
rotation to ij>l . In an equatorial orbit, the z coordinates are always 
identically zero in the integration. However, in an inclined orbit, numer- 
ical roundoff allowed $l to assume a slightly non-zero value, causing an 
erratic distribution of the angle between <j>i and after the first four 



output points . 
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In order to avoid the singularity in the angular representation, a 
pair of runs was done with non-zero initial values of .1, .2, and .3 de- 

grees for 6*, <f>l , and ipl in the rotating orbital reference frame. Runs at 
0° and 28° orbital inclination gave identical time histories of the angles 
measured in the orbital frami. . Figures 2.1.2, 2.1.3, and 2.1.4 give plots 
of the Euler angles 01, <f>l , and ipl for the 0° inclination case, and Eigures 
2.1.5, 2.1.6, and 2.1.7 give the angles at 28° inclination. Small differ- 
ence appear in the least significant digits. 

The last test run is the special case where the subsatellite precesses 
at the same rate as the orbital angular velocity so that the symmetry axis 
maintains a fixed orientation with respect to the orbital reference frame. 
The initial conditions are with equal to the orbital angular velocity, 
and a 28° orbital inclination. The angles 6> and <j>> remained constant 
during the run just as in the equatorial case run with the original version 
of the program. Eigures 2.1.8, 2.1.9, and 2.1.10 show plots of the angles 


vs. time. 
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Figure 2.1.2 0' [radians) vs. time [sec) in an equatorial orbit with 

initial conditions 6 0 « - .1°, * 0 ' = .2°, and %' = .3°. 
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Figure 2.1.5 0’ (radians) vs. time (sec) in a 28° inclination orbit 

with initial conditions 6 0 * = .1°, ^ 0 ' = .2°, and ^ • : 
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Figure 2.1.7 i|i’ (radians) vs. time (sec) in a 2b° inclination orbit with 
initial conditions 0 O » = .1°, $ • = .2°, and tl» 0 * « .3°. 
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Figure 2.1.8 0' (radians) vs. time in a 28° inclination orbit for the 

special test case where the precession rate is equal to 
the orbital angular velocity. 
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Figure 2.1.10 (radians) vs. time in a 28° inclination orbit for the 
special test case where the precession rate is equal to 
the orbital angular velocity. 
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2.2 Slack Tether Studies 
2.2.1 Introductory Remarks - 

Several advances have been made in the consideration of the slack 
tether (tether break) problem during the reporting period. 

Perhaps most significantly, the SLACK2 program described in the last 
quarterly report has been extended to SLACK3, a fully three-dimensional 
treatment : 

- The motion of the masses is three dimensional, 

- The deployment boom may point in any direction, out of the or- 
bital plane as well as in-plane. 

- The tether may be deployed in any direction (though it is still 
restricted to a straight line deployment) . 

- The Shuttle is allowed to rotate about an arbitrary axis at a 
specified rotation rate and/or rotational acceleration. 

- The vibration of the boom on release takes account of the various 
directions involved, e.g. the direction of initial tension due 
to the tether . 

SLACK 2 previously treated damping in the tether (when the masses 
caused a section to come into tension and "bounce") in an ad hoc fashion. 
A more correct treatment was developed and applied in SLACK3 . 

RECOIL, a high resolution program to model the initial loss of tension 
process following a tether break, was written. The results (for a per- 
fectly elastic tether) verify the initial conditions we have been using in 
SLACK2 and SLACK3 : the tether is recoiling at a constant velocity along 

its length, with no deformation. These initial conditions had been used 


Jh -c.,? 


on 
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the basis of some low resolution SKYHOOK runs and plausibility arguments. 

SAO, with Prof. Robert Hohlfeld of Boston University, has also begun 
consideration of slack tether dynamics, and the tether loss-of-tension 
problem, from an analytical viewpoint. These investigations are expected 
to provide valuable results not readily obtainable from computer modeling 
(e.g., the typical scale size of loop formation in the cut tether), to 
provide a check on computer results, and to guide the development of model- 
ing algorithms for partially slack situations. Prof. Hohlfeld has contrib- 
uted a section, 2.2.5, detailing the results so far. 


2.2.2 Treatment Of Damping - 

The algorithm for treating damping in SLACK2 was simple and ad hoc : 
When two of the masses in the model separated by the natural length of the 
(assumed massless) tether segment between them, they undergo a "bounce," 
reversing their momentum components along the segment joining them. Allow- 
ance was made for damping by reducing the velocity of the rebound by a 
constant percentage. This much is reasonable, simply from consideration of 
the linearity of the processes (see below) . What is ad hoc about our 
treatment is the way in which we determine what the percentage velocity 
reduction is. In the absence of any better knowledge, we simply input a 
percent value, which applied to all segments regardless of length. 

We have now put the damping considerations on a more physical level, 
so that we not only have the percentage velocity reductions for the various 
segments to proper scale, but compute them from an input constant charac- 


teristic of the tether material. 
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Although our results are certainly not new, we have not found them in 
a search of readily available sources in any form easily applicable to the 
situation at hand. Hence, we shall give a simple derivation from first 
principles. We wish to derive the dependence of < ne damping constant for a 
piece of tether material on the length of the piece, under some very gen- 
eral assumptions. Specifically, we assume that the material is viscous but 
not plastic . That is, there will be resistance to motion, but no permanent 
change in the equilibrium state. We assume a simple Hooke's law form for 
the elastic portion of the stress, and the viscous damping will add a term 
proportional to velocity. 

Consider a simple physical system: a length of tether fastened to a 
wall at one end and a mass M at the other. Neglect the tether mass. Let x 
be the extension of the segment past natural length. Then the equation of 
motion is 


Mx = - Kx - Bx 


( 1 ) 


where K is the Hooka's law spring constant and B is some damping constant. 
Now suppose we cut the tether in its center and place a small mass m 
between the pieces; denote the extension of the segment: between the wall 

and mass m by y. Then the equations of motion are twofold: 

" 

Mx = -k(x-y) -b(x-y) 

I my = — ky -by (2) 

+ k (x — y) + b (x - y) 

If we let ( = (x-y) , a bit of manipulation arrives at the equation 


= — tjkx 


%bx 


- ifa(x - 0 




&s •«C„ .3 V - • 


Mx 


O) 


Page 42 


Let m -+ 0 and compare the result to equation (1). If the accelera- 
tions in the third term of (3) remain finite, which we shall simply assume 
since this is not meant to be a rigorous proof, agreer. it between these two 
results for the motion of mass M requires that k = 2K, b = 2B. Generaliz- 
ing the argument , we have that 

k * J b « \ (4) 

where t is the (natural) length of the tether segment. For the spring 
constant k, the proportionality constant is simply EA, Young's modulus 
times the cross sectional area. For the damping constant b, we shall 
follow a remark in Bodley and Park (1983; p, 54) and denote the propor- 
tionality constant by C v . It seems reasonable that the energy dissipation, 
hence the resistive force, will be proportional also to the cross sectional 
area of the tether; hence, we define c v = Cy/A. The two material proper- 
ties E and c v are thus input to a program such as SLACK3 , which can then 
compute the spring motion parameters. 

We now need to determine how the velocity after a "bounce" in SLACKS 
depends on the relative velocity of the pair of masses before the "bounce." 
Consider only the velocity component along the tether; the orthogonal 
component cannot be changed by tether tension forces . Suppose we have 
equal masses m at either end of the tether (this won't be strictly true 
since SLACK3 allows the segment length to vary, and the masses depend on 
segment length), and a segment of length £. Let the separation of the 
masses be x+£, so that x is the stretch of the segment. The displacement 
of each mass will be hal f x, so the equation of motion (acceleration = 
force) becomes 


ijmx — 


— kx — bx 


(S) 
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Make the usual trial substitution x(t) — a 0 *, and write p — b/m , u> 2 = 
2k/m. Then the exponential will be a solution if 

a 3 + 2pa + w 3 = 0 (6) 

There are two linearly independent solutions corresponding to the roots a = 

~P ± \/W - w 3 of this characteristic equation. We consider here only a 

small damping a pproximation . i.e. p « w. The complications become more 
complex for larger damping, and the implications of extreme cases for the 
model are still being investigated; the results should be good enough for 
most purposes when P < w. The two roots then form a complex conjugate 
pair and the solutions are more familiarly written as exponentially damped 
sine and cosine terms. If the bounce begins at t=0, we have Initial condi- 
tions x(0) = 0, x (0) = v in , say. The cosine term will drop out and the 

solution becomes 

x (t) = e ^ cos ( w t) (7) 

where we have replaced %/ p 2 -w 3 by o>. The tether segment comes out of the 
bounce (i.e., loses tension again) when x (t) is next zero, i.e. at t — 
ir/w. Taking the derivative of (7) at this time we have 

v ou t = v in e'*/*' (8) 

The computation of segment damping parameter, and this small damping 
approximation for the velocity loss on bounce, have been implemented in the 
version of SLACK3 described in Section 2. 2. 3. 2 below. 

References to Section 2.2.2 

Bodley, C. S., and Park, A. C. , Analysis Tsthersd Satellite System Qn- 
bita_I Dynamic s for Selected Mission Profiles. Report TSS-83-ACP-065 , 
Martin Marietta Denver Aerospace, June 1983. 
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2.2.3 Full Three Dimensional Treatment - 

SLACK 2 restricted the motion of the masses to the orbital plane. The 
same physical model has been generalized to a three dimensional program, 
SLACK3. This involved the following modifications: 

- The equations of motion of the masses now allow for out-of-plane 
motion . 

- The deployment boom may point in any direction, out of the or- 
bital plane as well as in-plans. 

- The tether may be deployed in any direction (though it is still 
restricted to a straight line deployment) . 

- The Shuttle is allowed to rotate about an arbitrary axis at a 
specified rotation rate and/or rotational acceleration. 

- The vibration of the boom on release takes account of the various 
directions involved, e.g. the direction of initial tension due 
to the tether . 

- Randomization of tether segment direction is in a cone about the 
nominal tether direction, rather than just in the orbital plane. 

Some additional modifications have been made (or are being coded) , such as 
the improved damping model discussed in Section 2.2.2, and streamlining of 
the input . 

There are still some restrictions from a fully adequate three dimen- 
sional treatment. The two most significant seem to be 

- The air drag (wind) iJ still in the orbical plane. 

- The deployment boom is assumed to be attached to the Shuttle at 
the Shuttle center of mass. 
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2. 2. 3.1 Preliminary Mathematics - 

The equations of motion of a free mass in the orbiting coordinate 
system (i.e., relative to the orbit center « the center of mass) are 

I x * ' = 3x + 2y' 

y' ' = - 2x' - D (1) 
z” = -z 

where the time has been scaled by the orbital angular velocity (1 to r = flt; 
()' means d()/dr; D is a scaled drag parameter, assumed constant; and the 
coordinate system is as shown in Figure 2.2.1. The third equation (for z) 
has solution z (r ) = z (0) Binr + z' (O)cosr, and addition of these solutions 
to the program is trivial . 

Much the more complicated part, both conceptually and in the actual 
programming, has been dealing with the three dimensional geometry and vari- 
ous angles involved. Although mathematically elementary, there is much 
room for confusion and failure to coordinate conventions, so we will spell 
out in some detail what we are doing. 

Figure 2.2.1 shows the orientation of our reference coordinate frame, 
the standard frame centered at the center of mass (Shuttle) and rotating so 
as to keep the x-axis vertical. Also shown are the orientations of the 
pitch, roll and yaw rotations discussed below. 

Figure 2.2.2 points up a distinction (in a two dimensional case for 
clarity) that it is extremely important to be clear about: that between 
physically rotating a physical object or vector by prescribed angles, and 
then asking for the coordinates of the new o b j.e. c t in the original . {or . 
reference^- coordinate system ; and rotating the coordinate system in pre- 
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scribed fashion and then seeking the coordinates in the new system of an 
unmoving object . In both cases, the "input data" are the coordinates of 
the original object in the original system. We shall call these two situa- 
tions "physical rotation" and "coordinate rotation" for brevity. Brief 
reflection, perhaps while pondering Figure 2.2.2, reveals that these are 
inverse situations: i.e. , a physical rotation by a given angle produces 
the same numbers as a coordinate rotation by the opposite (negative) angle. 
We shall consider the physical rotations and their generating matrices to 
be the fundamental concepts . 

The three fundamental units from which we shall construct all our 
rotations and angles are pitch, roll and yaw. These operations are always 
taken relative to a coordinate system fixed in the Shuttle (body coordinate 
system) . They are somewhat easier to describe if we imagine the Shuttle to 
be oriented as in Figure 2.2.1, with nose pointing along-orbit and the 
landing gear pointed toward the Earth. Pitch by a positive angle then 
brings the nose down; roll causes a roll about the Shuttle's longitudinal 
axis clockwise as viewed by the pilot; and yaw by a positive angle swings 
the Shuttle's nose to the left. Unless otherwise specified, 6 will refer 
to pitch, <p to roll, and x/t to yaw. 

Each of these elementary maneuvers can be represented by a 3x3 matrix. 
These matrices are used to transform column vectors, e.g. for pitch 
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X 


cos 0 

- sin 0 

0 


X 

y 

= 

sin 0 

cos 0 

0 


y 

z 

1 

0 

0 

1 


. 2 . 


or in more compact notation 

?i = P o r 0 


(2) 


(3) 


Here, the subscripts 0 and X refer to before and after the physical rota- 
tion, and the r = [x y z] T are the coordinates in a fixed reference frame, 
coincident with the body frame before the pitch maneuver, of some physical 
object or vector which rotates with the Shuttle. Note that r is not a 
physical vector, but a column matrix of components. If the system in which 
these components are measured is not obvious , we may use notation such as 

Tlixys 


The roll and yaw matrices are similarly: 


Rp 


Y* 


cos <p 0 sin <p 

0 10 
— sin <p 0 cos <p 

10 0 
0 cos - sin 

0 si n^ cos 


( 4 ) 


(5) 


These rotations and matrices are defined so that the body coordinate system 
coincides with the reference coordinate system when the Shuttle is flying 
in a "normal" attitude, i.e. nose forward and head up. There may be some 
discrepancy in sign from the conventions used in, e.g., aircraft opera- 
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tions . 

These matrices are examples of orthogonal matrices, i.e. those for 
which the transpose is the inverse. For these particular matrices, replac- 
ing the rotation angle (e.g. 0) by its negative results in the inverse 

matrix. From the comment made above about the inverse nature of physical 
• rotation and coordinate rotation we see that the matrices for transferring 

the representation of a constant vector into a system rotated with the 
Shuttle will be just the transpose of those written down above; a general 
result for orthogonal rotation matrices. 

It now becomes relatively simple to compute the rotation matrix corre- 
sponding to successive rotations, e.g. a sequence of pitch/roll/yaw maneu- 
vers . In each case we suppose the rotation is described in the body coor- 
dinate system. First, perform a rotation described by matrix A, then one 
described by matrix B. These will transform a physical vector ?o to ?i to 
r a . Initially, the body system is coincident with the reference xyz sys- 
tem; after the first rotation, it will be coincident with a system, say 
x'y'z'; and finally, with Now, almost by definition, the coordi- 

l. nates of the final vector in the final system are numerically identical to 

| those of the original vector in the original system. I.e,, 

I ?J I ic*y’r" = ?0 I xyc 

, and we may then write immediately r 2 |xy B - A {rj| x v«'> = A <B rj| x - y »g*> 

I 

; = AB Policy*. This final expression is just in the form required for *’he 

l' total rotation matrix; i.e. , the total rotation is described by the matrix 

| AB. In general, a series cl£ notations.,, first ^ then EL. than will 

I produce a rotation described fey the matrix • The order is important, 

I 

! 


1 
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since matrices do not commute. 

With this background, the required computations are almost obvious. 

Any direction, e.g. boom or tether, will be specified by two angles, first 
a pitch and then a roll which are pictured as rotating the vertical (x) 
axis into the desired direction. The full rotation matrix is easily com- 
puted as PR, and the unit vector in the appropriate direction is given by 
PR [1 0 0] T <= [cos0cosV9 , sinflcosy? , -sin^s ] T . These unit vectors allow 
us to create suitable initial conditions. 

Dealing with the boom vibration and Shuttle rotation is slightly more 
complicated. First, suppose we have computed the boom's position and ve- 
locity Bo and B 0 relative to an unrotated Shuttle, To describe the Shuttle 
rotation we have the angles 6 and <p of the rotation axis, which allow us 
to compute a transformation matrix PRj and the rotation angle ^(t) as a 
function of time. Our prescription for performing the desired rotation can 
be stated; First, physically rotate the Shuttle with PR so that the body 
{-axis (vertical in the body frame) lies along the axis of rotation; sec- 
ond, rotate the Shuttle by ^ about the body vertical axis in a yaw maneu- 
ver, with matrix Y; third, rotate the Shuttle in the inverse manner [PR] T 
so that the body vertical axis is coincident again with the reference frame 
vertical. From the theorem given above on combining rotations, this is 
equivalent to using a single total rotation matrix 

A(t) = [PR] Y 0(t) [PR]* (6) 

that is, the Shuttle position at time t is given by 


Bi = A(t) B 0 (t) 


(7) 
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The velocity is found simply as 

Bi - A(t) B 0 {t) + A(t) B 0 (t) (8) 

where 

A(t) = [PE] [PR]* tf(t) (9) 

We still need to know the unrotated boom vibration, i.e. relative to 
the Shuttle. We already in SLACK2 compute the frequency and amplitude of 
the boom vibration, as detailed in a previous report. There, to get the 
amplitude we used the tether tension before the break based on vertical 
deployment. If the tether pulled away from vertical by air drag, and is in 
equilibrium, the tension is the same as if it had been vertically deployed 
with only the gravity gradient force} if the tether was in motion or held 
off vertical by thrusters, the force depends on factors and scenarios not 
input to SLACK3 . Thus, we simply use the vertical deployment tension, with 
an optional "fudge factor" to adjust to specific cases, To treat the three 
dimensional case we use the two unit vectors along the boom and tether, eg 
and e T , which we have computed as above based on the input angles. The 
displacement of the boom tip due to the deployed tether will be along a 
direction perpendicular to the boom, in the plane defined by the boom and 
tether direction. First note that the vector (egXei) is orthogonal to this 
plane. Then a suitable vector is 

Xp = (esX a?) X eg = Sg — (eg • ©g) eg (10) 

which we then normalize to form a unit vector 


e p = Xp / |Sp| 


(ID 
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(In the exceptional case that the tether and boom are collinear, x P is zero 
and cannot be normalized, but there is no displacement and we simply use 0 
in place of 6 P .) The displacing force, used to compute the boom deflection, 
will then be simply 


Fp = ■{Tension} (Si * ep) . (12) 

The initial, pre-break, displacement of the boom tip is then 

35 = E ’ 3W s ’ ■ < 13 > 

where b is the boom length and El its stiffness. The motion of the boom 
tip, prior to overall rotation of the Shuttle, is then 

B 0 (t) = B. + (AB) cos [2»rft] (14) 

and its velocity 

B 0 (t) = - 2jrf (AB) sin[2jrft] (15) 

- Ig El 

where B. = be B is the boom equilibrium position and f = ci y— — — 3 is 

frequency of the first mode of a clamped beam; we take values specified by 

NASA, El = 1.3 x 10 B lb in a , b = 849 in, weight W = 106.2 lb; the value of 

gravity is g = 385.06 ln/sec 2 , and e* is 0.56 cycle/sec. Rotating (14) and 

(15) with (7) and (8) , we specify completely the motion of the boom tip. 
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One other application of these concepts is to a problem that at first 
sight appears trivial, the generation of vectors randomized about a given 
direction. We do this in randomizing the directions of the tether segments 
when we generate the initial conditions . We first compute (random) angles 
6 with respect to the reference direction and p azimuthally about the ref- 
erence direction (typically 6 might be uniform with 5° mean and p uniform 
on [0,2?r) ). Then we create a vector with these angles about, say, the x- 
axis in the reference coordinate frame: 


s 


cos 6 

sin 0 cos p 
sin 6 sin p 


(16) 


Finally, we physically rotate this vector using any rotation which brings 
the x-axis into the desired reference direction. If the reference direc- 
tion is e T , we need a unitary matrix satisfying: 

- ®t UU T = I ( 17) 







I f we partition U B [ Uj t u a I ua ] , then we see immediately that u a must 
be e T , and {t.ij. , u a , u j) - must form an orthonormal set. This may be easily 
done by choosing any arbitrary vector t not parallel to Uj. (we use unit 
vectors along the coordinate axes for simplicity; we will never have to 
try more than two) . Then form Ui.xt and normalize it to get u a ; and then 
take u 3 = UiXUj. We may new perform our rotation, which will bring the 
segment into an angle 8 with respect to the tether direction; the effect 
on the azimuthal direction is not obvious, but it will leave the uniform 
distribution unchanged. 

2. 2. 3. 2 Slack3 Implementation - 

The above considerations have been implemented in SLACK3. The im- 
proved damping of Section 2.2.2 is also included, as well as improvements 
in the randomization processes. We are in the process of coding stream- 
lined input procedures. in which there are default values of many 
parameters and a menu system allows changing only those of interest. Only 
a few test cases have been run (which verify correct deployment of boom and 
tether, and boom rotation); the familiar side and rear view configuration 
diagrams are shown in Figure 2.2.3 for a sample run which includes boom 
rotation. SLACK 3 , with minimal further effort, is now suitable for running 


2.2,3. Sample res 
a of configuration 
ive configurations 
is a side view, 
-Diane" fioure . th 
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some complex cases in the next reporting interval to determine the effect 
of deployment directions, rotation maneuvers, tether remnant length and 
other quantities of interest, 

2.2.4 High Resolution Loss-Of-Tension Model - 

The initial conditions used by SLACK 2 , that the tether was recoiling 
with constant velocity and no compression initially, were based on some old 
SKYHOOK simulations with only a few masses, -and on some plausibility argu- 
ments. It was decided to check these assumptions by means of a lumped mass 
model similar to SKYHOOK, but simplified to the extreme to allow a large 
number of masses to be included. 

The model, embodied in a program RECOIL, is capable of handling up to 
100 masses with reasonable efficiency. The restrictions made are prima- 
rily '. 

- The only forces included are the internal tether forces. Specif- 
ically, no gravity gradient, drag, or Coriolis forces are in- 
cluded . 

- All motion is one dimensional. 

- All segments are the same length. 

- Once a segment goes slack we ignore any tension that may develop 
if it again becomes stretched. 

The equations of motions are solved with a simple fourth order Runge-Kutta 
integrator, starting from the moment the tether is severed. Eor simplic- 
ity, the problem is scaled so that the only free parameter is the damping 


constant . 
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Figure 2.2.4. High resolution loss of tension; development in time. This 
is a view of a three dimensional surface whose height represents the (nega- 
tive of the) recoil velocity. This is a function of the time and the 
distance along the tether, in the "horizontal" plane: the tether/boom 
attachment is at the right, the cut end at the left, and successive config- 
urations recede into the distance. Fifty mass model with RECOIL; zero 
damping. 


Sample results are shown in Eigure 2.2.4, which shows a three dimen- 
sional plot of the velocity vs. time and distance along the tether. Note 
the propagation of a sharp wave front, behind which the velocity is approx- 
imately constant. The state of the tether at the moment it goes completely 
slack is shown in Figure 2.2.5. It is seen that these detailed calcula- 
tions confirm the initial conditions we have been using. 

Attempts to run RECOIL with non-zero damping lead to sharply oscillat- 
ing results. There is probably a coding error in the program or a numeri- 
cal instability. We shall investigate this in the next reporting period. 
The efficiency of the approach used suggests the possibility of adding 
extra forces (e.g. gravity gradient) as well as damping. 

It is interesting to see how this lumped mass model relates to the 
tether considered as an elastic, continuum. By letting the number of tether 
segments become infinite, and scaling the ball masses and segment spring 
constants appropriately, we can derive a partial differential equation, 
which is (as expected) a simple wave equation. Using a standard low order 
finite difference approximation we arrive at exactly the same equations as 
in the ball -and-spring model. E xcep t that the cut end boundary condition 
is not treated appropriately. There are several ways to treat such free 
end (zero stress) boundary conditions which maintain the order of approxi- 
mation; simply considering it as another ball with the same mass as the 
Interior balls in the lumped mass model, however, is equivalent to treating 
the boundary condition with a ower order a pproximation than the inte- 
rior points. It is a truism of numerical analysis that the order of the 
total solution is strongly influenced by the weakest link. (This sort of 


consideration is hidden by the lumped mass approach of first discretizing 
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the physical system conceptually before writing down any equations, which 
is why lumped mass models should always be examined carefully.) By properly 
treating the free end we might substantially improve the accuracy of the 
model (which can be easily checked in the zero damping case) . The possi- 
bility of improving accuracy by, for instance, using a finite element ap- 
proach, is also appealing. 


It is also possible to derive the velocity and strain upon loss-of- 
tension by theoretical arguments. This solution is implicit in the mate- 
rial of Section 2.2.5, but not actually derived there. Thus, it is of 
value to write down the argument explicitly. We assume there is no damping 
or other energy loss. The loss of tension will then occur via the propoga- 
tion of a constant amplitude tension wave with velocity c = \/E/p . Since 

the wave is propogating along a characteristic (in the space-time plane) 
with no dissipation, it will remain constant in form and amplitude. Thus, 
the velocity of recoil after it passes will be the same at all points. To 
be concrete, suppose the tether is fixed at x = 0, is of natural length L, 
and has been stretched by a small fraction 6 . Then, where the coordinate s 
describes the natural length of tether from the fixed end, and x its actual 
position, what we have described may be written: 


x(s,0) 
dx (s , t) 

at 


~ (1 + 6 ) s 



(evenly stretched Initially) 


t < 
t > 



where V is the recoil velocity which we shall determine below, 
immediately have: 


( 1 ) 

( 2 ) 


Then we 


x (s , t) 


(1 + 5) s 

(1 + 5) - vj^t - L ~ 5 j 


t < 



(3) 




< 


t > 


c 
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Upon evaluation at t = L/c , when the wave reaches the point of attachment, 
this is simply 

x(s,L/c) ** [l + * “ “ ] (4) 

The final step is to compute the recoil velocity. Our assumption of no 
dissipation allows us to simply equate the elastic energy stored in the 
tether segment with the kinetic energy after recoil. The spring constant 
of the tether as a whole is k = EA/L, and the force upon extension to L+u 
is just F - -k u. Integrating Fdu from 0 to u = 6L, we get the total 
elastic energy stored when the tether is stretched to length (l+tf)L, 
E«u§tic ” Jjk(£L) 2 ~ ijEAL6 2 . Equating this with the total recoil kinetic 
energy, Efciiwtic = (mass) (velocity ) 2 = ij{pAL)V a , gives 



and then evaluation of (4) gives simply x(s,L/c) = s. I.e., there Is. ns 
residual strain . Thus , the initial conditions we have been using in 
SLACK2/3 are exact under the assumption of zero dissipation. 


-* JSk- •• 
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2.2.5 Analytical Studies Of The Slack Tether Problem - 
2. 2. 5.1 Introduction - 

A program of analytic studies of tether dynamics pertinent to treat- 
ment of Lhe slack tether problem has been undertaken. While it is not 
realistic to expect that a complete understanding of the slack tether prob- 
lem may be obtained solely by analytic techniques, compelling reasons 
motivate analytic studies as an extremely useful adjunct to numerical mod- 
eling of the slack tether. Among these reasons are: 

- Analytic calculations may give insights into the physical proc- 
esses important for tb.e dynamics of the slack tether. Some of 
these processes may be obscured by the complexity of the numeri- 
cal simulations. 

- Numerical modeling of the slack tether using the present ball- 
and- spring techniques is computationally intensive. Analytic 
studies will help to direct the numerical modeling, and perhaps 
to suggest new, efficient modeling techniques. 

- Analytic calculations will aid in the validation of the results 
of the computer -generated numerical models. 

- Information generated by the analytic work will frequently com- 
plement the results of numerical simulations. 

The analytic calculations will be initially focussed on processes oc- 
curing just as the tether is going slack, i.e. going from a tensioned to an 
untensioned state. In this context, it is appropriate to treat the tether 
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as a one-dimensional system, as long as the "loops" one would expect to 
develop in the slack state do not grow to large amplitude in the transverse 
directions . 

The desired results of the initial theoretical effort is a calculation 
of the characteristic size for loops or wrinkles in the tether after it has 
gone slack. The characteristic size is of immediate importance for the 
safety of Shuttle operations in a slack tether situation. We expect that 
loops of small scale size (say on the order of 1 meter) will be much less 
apt to entangle the Shuttle than larger-sized loops. We shall outline 
below the progress being made toward this and other objectives. 

The first step in this program of calculations is to determine condi- 
tions pertaining to the development of patterns of tensioned and slack 
(untensioned) regions along the length of the tether. For simple cases, 
such as a tether break near the tethered satellite, it is obvious that the 
tether will go slack through the mechanism of propagating a longitudinal 
(negative amplitude) stress wave along its length. More complex situations 
may be envisioned as a superposition of (positive and negative) finite- 
amplitude longitudinal stress waves propagating along the tether. An imme- 
diate result (which will not be given in detail here) is that when the 
tether fails, it is possible in principle, to have a nearly arbitrary dis- 
tribution of tensioned and slack regions along its length. In particular, 
It is possible to construct initial conditions for the subsequent dynamics 
after tether failure consisting of a tensioned region bounded by two slack 
regions . 

Such a discussion of longitudinal stress waves propagating through the 
tether suggests we apply the method of characteristics to the dynamics of 
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the tether as it is going slack and for times shortly thereafter (Morse and 
Eeshback 1953, Courant and Hilbert 1962) . Recall that a characteristic of 
a t«ave equation (hyperbolic partial differential equation) is the locus of 
points in time for the propagation of an Infinitesimal amplitude wavelet. 

A cable or tether going slack can be modeled (for times before the 
formation of well-developed loops) as a column capable of supporting a 
great deal of force in tension, and of supporting a comparatively small 
(but nonzero) force in compression (Landau and Lifshitz 1959, §§17, 20). 
The formation of transverse bends in the tether which will subsequently 
develop into loops is a generalization of the well -studied problem of col- 
umn buckling. 

Since the compressive force the tether can support is very small, a 
number of effects will have to be considered which are usually neglected 
for models of the tether in tension. While the tether in tension can be 
modeled as a string with a Hooke's law constitutive (stress-strain) rela- 
tion, that slack tether will require a nonlinear viscoelastic or viscoplas- 
tic constitutive relation describing its resistance to transverse bending. 
Also, mathematical treatments of the slack tether will require a careful 
treatment of the tether inertia. 

Recall that the dominant "mode" of column buckling for a column sup- 
porting a large compressive force and with negligible Inertia is an approx- 
imately half-sinusoidal displacement over the entire column. Qualita- 
tively, in this instance of generalized buckling, the tether inertia will 
favor wrinkling or buckling on smaller scales and the resistance to bending 
will tend to favor larger length scales. Consequently, we can hope to 
identify a fastest-growing scale size for buckling which will characterize 




loop formation in the slack tether, 
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2. 2. 5. 2 Motion Of A Medium Under Its Own Elastic Forces - 

We shall develop a general equation of motion for a medium under the 
influence of its own elastic forces. A special application will be made to 
longitudinal stress waves in £ cable, such as the tether of the TSS , in the 
next section. While the general formalism is not necessary for this first 

simple application, it is valuable to introduce this notation for later 

calculations pertaining to the slack tether. 

Following Morse and Feshbach (1953) , define the stress dyadic 

I = E x i + F y j t F B k a) 

where i, j, and k are the usual unit vectors in rectangular coordinates. 

The force across surface element dA is then S • dA , 

The stress dyadic is related to the strain dyadic defined in (3) by 

I = Al 6ll + 2/i6 (2) 

where I is the unit dyadic. Equation (2) is immediately apparent by trans- 
formation to the principal axes. The quantity, /*, is the shear modulus of 
the solid, ^ (3A+2^i) is the bulk modulus of an elastic Isotropic medium,, 
and /i (3A+2/i) / (A+fi) is the Yeung’s modulus. The strain dyadic may be writ- 
ten in terms of the displacement, s, of the solid from its equilibrium 


position , 
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6 


1 

2 


(V§ + sV) 


(3) 


The equation of motion for the solid then, under the influence of its own 
elastic forces is 


P 


9 a s 

dt a 


V* [AIV* § + /iVs 


^sV] 


= (A+/z)V(V-s) t /iV- (Vi) 


= (A+2/-J) V(V>s) - ^iVx(Vxs) 


(4) 


The form of the equation of motion suggests that at least part of the 
vector displacement, s, may be expressed in terms of the gradient of a 
potential, rp. Let 


s = V\p 


and so 



dhp 

dt* 


(5) 


where 


C c 2 = 


A+2 fi 
P 


which is the usual wave equation, a hyperbolic partial differential equa- 


tion . 
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If the shear modulus is zero, the wave speed squared is Just the bulk 
modulus divided by the density, as expected. The longitudinal character of 
these waves is immediately apparent, since if the solution of equation (4) 
is taken to be a scalar, the dyadic operator 

D = V(fyi) 


is symmetric, and the rotation dyadic 

t « - |{Vxs)xl = 0 

and so D = 0, the strain dyadic. For such a solution, there is no twisting 
of the medium, only stretching and squeezing, and so the wave is purely 
longitudinal (compressional) . Therefore, for purposes of this discussion, 
it is possible to treat longitudinal waves in a Cable as a special case of 
longitudinal waves in general solid media without loss of generality. 


2. 2. 5. 3 Longitudinal Stress Waves For Simple Viscoelastic Constitutive Re- 
lati. s - 

The general case of a time-varying tension in a cable or tether can be 
treated using the method of characteristics as the superposition of a num- 
ber o finite-amplitude longitudinal stress waves {Morse and Feshbach, 
1953; Courant and Hilbert, 1962) . For the case of a completely linear 
elastic constitutive relation (Hooke's law) , wave propagation is ncn-dis- 
persive. However, for the case of loss of tension in the tether, visco- 
elasticity must be taken into account. 
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We shall examine wave propagation in the context of two models of 
visco-elastic body strain, the Maxwell model and the Kelvin model (Drucker 
1967) . More complex models may perhaps b' ~equired ultimately to account 
for details of slack tether dynamics, but ter the present we shall use 
these models to illustrate changes in elastic wave propagation by small 
perturbations of viscosity. 

In the Maxwell model , the strain is taken to be the sum of elastic and 
viscous components, i.e., 


e = + e v 


and the elastic and viscous stresses are constant at all times, 


(7 V - <x a ~ a 


The elastic stress and strain are related by 


e° = <t/E m 


and for the viscous component, 


de v 

dt 


= o/C M 


independent of all prior history. The equivalent mechanical system for the 
Maxwell model is a spring and dashpot connected in series. 
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We can write a differential equation for the total strain, e, 


de_ _ de^ de^ _1_ d£ jt_ 

dt dt + dt E m dt + Cm 


(Note that treatment of the strain, e, in this section is equivalent to the 
treatment of the displacement, s, of the previous section.) Since there is 
no restoring force associated with the displacement of the "dashpot," i.e, 
the viscous strain, the perturbed wave equation for the Maxwell model is 


d i 

3x 2 


6 ” = 




+ e v ) 


but d« v /dt = a/ Cm, and so the equation can be rewritten as a wave equation 
for the stress in the cable, 


d 2 a _ _1_ d 3 _ B h da 

C c 2 ' dt 2 ' T ~ C c 2 C„ 3t 


We shall consider the viscous term as a small perturbation on the 
hyperbolic partial differential equation and attempt to find a solution of 
the form 


a — f(x±c c t) exp ^-a|x-x 0 |j 

a wave spatially damped from some initial point of propagation, x„. Sub- 
stitution of this form into the wave equation gives, to lowest order in 
viscous terms , 

Em 

2C e C M 


a = 


x - ... . v i- .-r » v ' ’ * -^-v.-,:. ».-■■ ^^-■■'rr'T^Tx.trrir 
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which relates the model parameters to the damping parameter of the wave. 
Therefore, treating a damped wave in the context of the Maxwell poses no 
difficulties or inconsistencies. 

A similar treatment may be undertaken for the Kelvin viscoelastic 
model. In the Kelvin model, the stress is at all times the sum of the 
stresses in an elastic and in a vinccuB element, 

a — (T m + a v 


and the strains are equal , 


£• = e v = € 


with 


ff* = E u c 


and 


<r v = C k 


de 

dt 


Independent of the history of the system. The equivalent mechanical system 
of the Kelvin system is a spring and a dashpot in parallel. The wave 
equation for the Kelvin system may, after some minor manipulation, be cast 
in the form of a wave equation for the strain, 


(JL 

1 d* \ 


ji 

(IL . 

_ 1 _ d±\ / 

\dx z 

C c 2 dt z) 

£ - - 

e k 

dt 

1 

Vdx2 

c c 2 at *) 


which manifestly describes the propagation of a damping wave with damping 
coefficient , 

Ek 


a ~ 


C c C K 


. • 
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Generalization to more , . mplex viscoelastic relations is straightfor- 
ward. It can be seen on the basis of these calculations that we can treat 
the propagation of tensioned and slack regions in a viscoelastic tether in 
a tractable fashion as a perturbative extension of the method of character- 
istics . 



2. 2.5.4 Analytical Studies Status - 

We have begun a series of calculations relevant to the slack tether 
pr-oblem when tension has been lost. At present, this work has been con- 
cerned with the situation in which the tether is going from a tensioned to 
a slack state. Further research will address dynamical problems of the 
slack tether. 

This program of calculations has begun from the application of the 
method of characteristics to the dynamics of longitudinal stress waves for 
a tether with a purely elastic constitutive relation. This work has first 
of all verified analytically the approximate validity of the constant ve- 
locity initial condition used in simulations of the slack tether in the 
program SLACK2 . Further, It has been demonstrated that for a tethe - with 
an elastic constitutive relation, arbitrary patterns of slack and tensioned 
tether segments can be built up by superposition of longitudinal stress 
waves. In particular, the important result has been obtained that is pos- 
sible to have a bounded, tensioned tether segment bounded by two slack 
segments, the boundaries moving with the velocities of characteristics for 
the hyperbolic equation for longitudinal stress waves. 
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Work is underway developing the dynamics of a slackening tether exhib- 
iting a viscoelastic constitutive relation, e.g. a Maxwell or Kelvin rela- 
tion. As the tension of the tether diminishes, eventually viscous effects 
must dominate the dynamics of the tether. This work has shown, as a per- 
turbation upon the hyperbolic equation for the elastic constitutive rela- 
tion, that the stress waves do not exhibit a nonlinear steepening for phys- 
ically reasonable constitutive relations. 

The work outlined above is of Importance for defining the initial 
conditions for a semi-analytic treatment of the slack tether problem. As 
outlined by Landau and Lifshitz (1959) , a slack tether can be modeled as a 
rod with small (not necessarily positive) tension with inertial forces and 
resistance to transverse bending. A taut string is the opposite limit in 
which a positive tension dominates the dynamics. The advantages of this 
type of mathematical treatment are several: 

1) A relatively consistent formalism may be developed applicable to 
both tensioned and slack tether segments. 

2) The beginnings of loop formation may be treated in a manner similar 
to the problem of a column buckling under compression. This would allow 
identification of the scale size of loop formation in the slack tether. 

3) No nonphysical mathematical singularities are introduced into the 
behavior of the tether at zero tension. 

4) Provided deviations from a linearly extended tether remain small, 
tether dynamics will be accessible by relatively standard mathematical 
techniques . 


~«CV.3v 
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2.2.6 Concluding Remarks - 

During the reporting ; 'od SAO has substantially enhanced the simula- 
tion program SLACK2/3 : The treatment of viscoelastic damping has been 
improved, and the program has been extended to full three dimensional capa- 
bility. In the next reporting period we expect to finalize the coding of 
SLACK3 (simplifiad input procedure and a few other minor itnms) and use the 
program to conduct a series of operational mission simulations relevant to 
safety issues. We shall also investigate the scale length for critical 
damping in the tether material and the implications for the bounce calcula- 
tions in the SLACKS model . 

An efficient, high resolution (yet strongly restricted) program for a 
partially slack ball-and-spring model has been developed. The program has 
confirmed the Initial conditions being used by SLACK2/3 (which assumes the 
tether has become completely slack and hence cannot follow the initial 
loss-of-tension) . These initial conditions have also been confirmed ana- 
lytically in the undamped case. In the next period we hope to extend the 
program to include damping. Other possible extensions include modeling 


with the gravity gradient force, and 2 or 3 dimensional modeling with Cori- 
olis and drag forces; and proper treatment of the free end boundary condi- 
tion from numerical analysis of the governing partial differential equa- 
tion. Perturbation (stability) studies might also be feasible. 

Analytical studies of the slack tether have been initiated, with the 
intent both of supporting the SLACK2/3 simulations and producing signifi- 
cant results not amenable to direct modeling. In the next period, we 
expect these studies to provide a good idea of the scale of the wrinkles or 
loops resulting from the initial buckling as a cut tether goes slack. This 
result will be relevant to safety studies and will also suggest how initial 
conditions for SLACK2/3 should properly depart from purely symmetric re- 
coil, a consideration now handled by ad hoc randomization. 

A continuing frustration, one which Is becoming critical for soma of 
the projected effort, is the lack of hard data on tether material proper- 
ties beyond the simple axial stiffness (EA) . We have asked Martin Marietta 
Corporation to send us what data they have on tether properties , and we 
also intend to request samples of tether from MMC, upon which we may per- 
form some s_ pie experiments . 
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2.3 Development Of A Computer Code To Model The Electric Potential Around 

A Severed Tether Immersed In A Plasma 

2.3.1 Introductory Remarks - 

In the Quarterly Report ttl, we showed that the prolate spheroidal 
model of the broken tether was very inaccurate at distances close to the 
wire. We are now developing a numerical model which represents the wire as 
a small cylinder. The program uses the Poisson’s equation to calculate the 
electric scalar potential on an axially symmetric three-dimensional grid. 
The grid lines are closely spaces’ near the end of the wire but the spacing 
increases geometrically as one moves away from the end. The lines parallel 
to the wire are spaced <3.1 mm apart inside the wire and within one wire 
radius. Then the spacing is increased geometrically as the distance from 
the wire is increased. Eor example by choosing a geometric factor of 1.2 
(each grid spacing is 20% larger than the preceedlng one) and by adopting 
for the first grid spacing a value of 0.1 mm, then 100 grid lines will span 
a distance of 41.4 kilometers. The radial grid lines are, likewise, 
closely spaced near the end of the wire, and their spacing increases as the 
distance from the end increases. 

While this approach could, In principle, be used to model the entire 
tether, we decided to develop and use a more efficient model. Rather than 
numerically modeling the entire length of the tether, the voltage at the 
grid points within a region spanning a space of 10 to 100 meters from the 
end of the wire is initialized by a prolate spheroidal model for the wire 
immersed in an otherwise uniform electric field created by the motion of 
the wire through a uniform transverse magnetic field (VxB) . The plasma 
density is initialized at this potential and then the trajectories of a 
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number of "computer particles" each of which represents many real particles 
is calculated. At each time interval the electric charge density on the 
grid is calculated and a discrete implementation of Poisson's equation is 
used to numerically estimate an updated electron potential. In this man- 
ner, the motions of the electrons and ions near the end of the wire can be 
modeled. While we have made considerable progress in developing this com- 
puter program the code is not yet operational . 

2.3.2 Generating The Coordinate System By Conformal Transformation - 

Prelate Spheroidal Coordinates are the most useful coordinates for 
calculating the electric field around a long wire. Moon and Spencer give a 
most convenient transformation for generating an orthogonal set of hyper- 
bolas and ellipses; 


z — c cosh{w) 


( 1 ) 


where 


V.. 
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z =« x+iy 

(2) 

w = rj + 10 

(3) 


Given a planar grid in the w = tf*l 0 coordinate system this transformation 
maps the straight lines »j = constant and 6 <= constant onto ellipses and 
hyperbolas in the z = xtiy coordinate system as shown in Figure 2.3.1. 
Rotation about the major axis of the ellipses generates the prolate sphe- 
roidal coordinate system. fhe purpose of this Section is to derive the 
relationships required to perform calculations in this interesting coordi- 
nate system. 







2.3,3 Elliptical Coordinates 


Substituting (2) and (3) into (1) we find that 


x*-iy = c cosh(r/ + i 0) 


= - + e'”' 10 ) 

2 


“ e'Hcos 0 + isin0) + e -, J(cos0 - isintf) 


= S. ( e r» + e -'j) C os0 + i (e" - e'^sinS 

O 


from which wb see that 


x = c coshr? cos6 


y = c sinhfj sin0 
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Therefore , 


c coshf ) 


+ ( (■■■ ^ = cos a 0 + sin 3 0 — 1 

J yc sinh »/ J 


Thus we have the equation for an ellipse where the major semi-axis is given 


a = c coshfj 


and the minor semi -axis is given by 


b = c sinhr; 


(9) 


From the definition of sinht; and coshrj it follows that 


cosh 3 »j - sinh a »7 = 1 


( 10 ) 



and so 


a a - b 3 = c 3 ( 11 ) 

At this point it is useful to introduce two vectors which extend from the 
two foci at x - ±c to any point on the x-y plane; let their magnitudes be 

ri = y/\x-c) 2 + y2 (12) 


and 


rj = \/ (x+c) 2 + y 2 


(13) 




Recalling that 


it follows that 




x = c coshtj cos0 


y = c sinhfj sin# 


( 5 ) 

( 6 ) 


■ S ~ ■ ’V'*' '■VTBr® 



Page 82 


ri = /(c coshij cos 0 -c) a + c 8 sinh a f? sin 2 0 


ri = cV , cosh a »j cos a 0 - 2 ccshtj cos 0 + 1 + sinh 2 »j sln a ( 


r x = cVcosh *ri cos *6 - 2 coshi? cos 0 + 1 + sinh 2 rj (l-cos 2 0 ) 


ri = c y/coe^d (cosh a t> - sinh a » 7 ) - 2 coshff cosJ + cosh 2 !? 


r-„ = c y/cos a 6 - 2 coshtj cos 0 + cosh 2 i 


rx = cy/~(cos$ - cosh 77 ) 2 


rx = c| cos 0 - coshtj 


But cosh*j > 1 5 c.osd , therefore 


rx = c (coshij - costf) 


(14) 


Comparing the definition of rx and r* one sees that the same logic used to 
evaluate r x in terms of r) and 6 results in 


r 2 = c (cosh ?7 + costf) 


(15) 


Therefore 


rx + r a = 2 c coshtj 


( 16 ) 
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rj " 1*1 “ 2c coed 


(17) 


Two new parameters are now Introduced which simplify the equations 


e = cosh, = 


(18) 


and 


ft = cos’? = ra ri (19) 

It was shown earlier that constant £ = cosh, defines an ellipse. Now we 
shall see that constant ft — cob$ defines a hyperbola 

[rsb?] ' [rfcj - c0 - fc - - 8lnh ’’' ■ 1 < 20 > 

The family of hyperbolas is orthogonal to the family of ellipses. This 
results fron the conformal transformation whi..h translates squares in the 
(,,d) coordinate system into curvilinear squares in the x-y coordinate sys- 
tem. This is easily verified by testing and mapping function to determine 
if the Cauchy-Riemann conditions 

§* - §y 

$Tf 39 ’ 36 " 3, w 

are met. Again recalling that 

S x = c cosh, cosd 
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y = c sinh?? sin0 


it follows that 


Bx . . 

•5— = c sinhfj coed 
dt] 

Bx 

— = -c cosh?? 

Bind 

(22) 

^ = c sinh?7 cob0 

= c cosh?? 
07] 

sin# 

(23) 


Thus we sea that the Cauchy-Rieraartn conditions are met and the set of 
hyperbolas are orthogonal to the set of ellipses as illustrated in Figure 
2.3.1. Rotating this coordinate system about the major axis of the ellip- 
ses generates the prolate spheroidal coordinate system j.n 3 dimensions 
H, and rp , It is conventional to rotate the ellipse with its major axis 
coincident with the z axis of the cartesian coordinate system. To accom- 
plish this one must replace x by z and y by R = \/x * + y 2 in the foregoing 
equations . The equations 


x = c cosh?? 

coH0 

(24) 

y — c sinh?? 

sLnO 

(25) 


are replaced by 


z = c cosh?? cob8 


(26) 


R = c sink*} sin# 


(27) 
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x “ R cos ^ 

(28) 

y = R slnV> 

(29) 

Therefore the relationship between prolate spheroidal 

coordinates and car- 

tesian coordinates is given by 


x = c sinh, sin# cost p 

130) 

y = c sinh, sin0 sinV> 

(31) 

z = c cosh, cos 8 

(32) 

or 


x = cy/£ s -l cost j> 

(33) 

y = cy/t 2 -! y/l-fi 2 sLnip 

(34) 

z “ c £ t* 

(35) 

where 


? = cosh, - ri 2c ra 


= cos * ~ 2c 
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ri = %/x a + y 2 + (z-c) 2 


r 3 r - \/x a + y 2 + (z+c) 2 


2.3.4 The Gradient, Divergence, And Laplacian In Prolate Spheroidal Coor- 
dinates - 

2. 3. 4.1 The Metric Coefficients - 

To perform field calculations in prolate spheroidal coordinates it is nec- 
essary to calculate the gradient, divergence and curl in that coordinate 
system. The easiest approach to this problem is to derive the gradient, 
divergence and curl in generalized coordinates and then to evaluate the 
metric coefficients hi, hj, and hj using the transformation from prolate 
spheroidal to cartesian coordinates derived in Section 2.3.3. 

To calculate vector relationships in generalized coordinates one must 
include the metric coefficients. In cartesian coordinates a differential 
length ds is given by 


ds = >/ "(dx) 2 + (dy) 2 ♦ (dz) 2 (1) 


In cylindrical coordinates ds is given by 

ds = y/ (dr) 2 + (rd0) 2 + (dz) 2 


not 


de = y/ (dr) - + (d0) 2 + (dz) 2 


WTO-7- 




>; * • fc* * * 
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In any orthogonal system, 


(ds) 3 = (h x dqi) 3 + (hjdqj) 3 + (h 3 dq 3 ) 2 


(2) 


where hi, hj and h 3 are the metric coefficients, 


Now to determine the metric coefficients we consider the following total dif- 
ferentials: 


dx 


+ feK + (^H 


dy + (^) dq2 + (^ j ) dq 3 

dz °(^ dqi + (^ dqz + (^) dq 3 


( 3 ) 


(A) 


( 5 ) 
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where q^, q£, and q^ are Lue generalized coordinates. 

Squaring these terms results in: 

<dx)2 " d ' n ] 2 + d ^] 2 + [ i ^ d< ^ 2 + cross terras (6) 

(dy) 2 = dt »i] 2 + [?q^ d( J 2 + ^i ] 2 + cross terms <7) 

(dz) 2 „ j^L dq J 2 + j^- d q 2 J 2 + ^ dq 3 ] 2 + cross terms (8) 

In any orthogonal coordinate system the Pythagorean theorem demands that the 
cross terras must vanish. Thus, assuming that our generalized coordinate 
system is orthogonal, it follows that 

"'■[(*9 '•($)■* Us) '>-.>■ 

Iffi) ■* ft ) 2 * ($) 


( 9 ) 
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Comparing equation (9) to (2) we 

see that: 


«; - m • * K) ■ 

+ (^) 2 ■ 

(10) 

Now to evaluate the metric coefficients for prolate spheroidal coordinates we 
make use of equation (10) and the coordinate defining equations obtained by 
the conformal transformation defined in Section 2.3.3. 

x = c sinhn sin0 cosij; 

o 

I A 

a 

A 

8 

(ID 

y ° c sinhn sin0 sin^ 

o i 8 i n 

(12) 

z => c coshn cos6 

o <. ^ £. 2tr 

(13) 


where 


n. q, c 0. q-, - 


The metric coefficients are calculated as follows : 

-Hff) 2+ (i*) 2 + (l0 2 

2 2 2 

n (c coshn sin0 cosif>) + (c coshn sin0 sin^) + (c sinhn cos0) 

2 2 2 2 2 2 

a c cosh n sin 0 + c sinh n cos 0 

D (1 + sinh^n) sin^B + c^ sinh^n cos^0 

= c /(sin 2 0 + sinh 2 n) (14) 

A - (fr) 2 + (Sf) 2 + (If) 2 

2 2 2 

= (c sinhn cosQ cosifO + (c sinhrj cos0 sirt/O + (c coshri sin6) 

2 2 2 2 2 2 

a c sinh n cos 9 + c cosh n sin 0 
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h 2 . hj 


h 2 . [2*\ 2 + [Sz\ 2 + /2 e\ 

n 3 \ty) V^/ V^M 


2 2 
= (c sinhn sin0 sinty) + (c sinhn sin8 costj)) 


=> c sinhn sin0 


2, 3, 4. 2 The Gradient, Divergence and Laplacian in Generalized Coordinates 
The gradient in generalized coordinates is given by: 

h l 8q l h 2 3q 2 h 3 9q 3 


where u^, u 2 » u^ are the unit vectors of the generalized coordinate system. 

The divergence of a vector is defined as the limit of the total flux emerging 
normally from a closed surface divided by the enclosed volume as the volume 
approaches zero. Thus 


_1_ f A“ < 
AV JAV 


where 


A = A^u^ + A^U2 + A^u^ 


Consider a differential curvilinear cube centered at (q^, q 2 , q^) as shown in 
Figure 2.3.2, The area of a differential square in the q 2> q 3 plane at the 
center of the cube is l^h^Aq^q^. The area of the end of the cube located at 
q^ + Aq^ is therefore given by Taylor’s theorem as 


larCTi 


„Kz 


Jl 
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Aa 


r g I 

Aq 2 Aq 3 ^ (h 2 h 3 ) — ) 


It follows that the normal flux through this surface Is given by 
A$ ■= A^'AS 

° At J 2 Aq 3 ^ h l h 2 A l + 9q^ ^ h 2 h 3 A l^ ~2~] 


Adding to this the flux through the opposite side of the cube in the -u, 

— ► 

direction, we see that the total flux exiting normally in the,u^ direction 
is given by 

A ''l " f h 2 h 3 + Sq^ < h 2 h 3 A l' T~ "] Aq 2 Aq 3 

- [ h 2 h 3 + 4^ (h 2 h 3 A l > i -r)] Aq 2 Aq 3 
But the volume of the curvilinear cube is 


AV » h 1 h 2 h 3 Aq 1 Aq 2 Aq 3 


Therefore the contribution to the divergence is given by 


(V.t) 

v 'l AV+O AV 


b [4; <h 2 h 3''l > 


(20) 


Applying this procedure to the other two axes and summing the results, we see 
that the divergence of A in generalized coordinates is given by 


v4 


tl l h 2 la 3 


[A 


(h 2 h 3 A^) 


9 

9q 2 


(hjh 3 A 2 ) 


9q 3 < h l h 2 A 3 } 




( 21 ) 


Now that we have the equation for the divergence and the gradient in generalized 
coordinates, it is clear that the Laplacian is given by 
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2. 3. 4. 3 The Gradient And Laplacian In Prolate Spheroidal Coordinates - 


Substituting (14), (15) and (16) into (17) gives the gradient operator 
in prolate spheroidal coordinates 


3* + . _ 9 . gi , (23) 


^ ** o'ij in *9 + sinh 2 n ^ c/sin ! 0 + sinh 2 n 99 c sin0 dip 


To calculate the coefficients of the Laplacian we utilize the following facts 


^1^2^3 ^ + sinh^n) slnhr) sin0 


(24) 


^2 h 3 

~r — = h, = c sinhn sin0 
h, 3 


(25) 


^1^3 

— r — = h_ ° c sinhn sin8 

H-n J 


(26) 


^1^2 _ c (sin^8 4- sinh^n) 
h, sinhn sin8 


(27) 


Introducing the prolate spheroidal metric coefficients into the Laplacian in 
generalized coordinates results in a cumbersome equation. It is clearer to 
deal with each term of (22) individually: 



- ^ (■= sinh i sin8 !£) 

c c sin8 | sinhn |“!r + coshn 



90 ( c 8inhr) sin9 30 ) 

c sinhn ^sin8 -|g^- + cos0 



(28) 


( 29 ) 
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JL ( hlhz ii.'l J_ ( (sin 2 e + sinh 2 n) 3£ 
3q 3 \ h 3 3q 3 / “ 


sinn sinQ 3i^ 


» (sln 2 9 + sinh 2 q) 3 Z (ft 
sinhn sin9 


Substituting (24) through (30) into Laplace's equation (22) results in: 


V 2 <+> 


c z (sin z 9 sinh z ri) sinhn sinG 


sin0 ^ 


3 z <t> 


sinhn rn T + coshn an 


an / 


+ 0 + cose ||) + sa£+ ga'-i $. ] 


which simplifies to: 


v 2 * = \ i!i + cothn iw + + r0 >.Q lil 

9 c*sin*e + sinh i n ) [ar^ + othr| 3n 39 2 + * 0W 39 j 


3 2 <f> 


c 2 sinh 2 r| sin z 0 3^ 


( 30 ) 


(31) 
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2.3.5 Solution Of Laplace'S Equation For The Electric Scalar Potential In 
Prolate Spheroidal Coordinates - 

It was shown in Section 2.3.4 that Laplace's equation in prolate 
spheroidal coordinates is given by 


v 2 4> 


c 2 (sin 2 0 + sinh 2 n) 


[0 + c ° thr ' In + 0 + cot0 ll] 


+ 1 _2!i 

c 2 sinh 2 p sin. 0 9 i|j 2 


( 1 ). 


We solve this equation by the technique of separation, of variables. Assume 
a solution of the form 


<j> = H(n)0(em^) 


(2) 


Substituting (2) into (1) and multiplying by c 2 sinh 2 r) sin 2 0/H04' results in 


sinh 2 n sin 2 


sin 2 0+sinh 2 


£ ( l d 2 H 
n \ H W 


L dH , 1 d 2 0 , - 

+ cothn ^ + qmt + cot0 


d0 1 , 1 d 2 f 
d0 ) V dip 


0 (3) 


The last term is a function only of l|* which does not appear in the rest of the 

1 d 2 ^ 

equation. In order for the equation to hold for any value of ip, — must be 


¥ dlj/ 2 


constant. By convention we choose this constant as -m , therefore 


d 2 ¥ 


di[i 


T + m 2 ^ = 0 


w 


From the geometry of the problem we know that (4) must have a single valued 
trigonometric solution; therefore, we deduce that 

=» A sin mtp + B cos mi|/ (5) 

and that m is a non negative integer. 


C-3- 




js* 


Equation (3) may now be rewritten as 


1 (d 2 H . .. dH \ 1 (d 2 G 

H {df? + cothn dn / + 0 IdP" 


+ cot9 


d©\ _ m z (sin z 9 + sinh 2 p) n 
d9 / sinh^n sin 2 9 = 


( 6 ) 


Now the rightmost term is easily shown to be 


m 


^sinh 2 n + sin^l 


(7) 


Therefore equation (6) is easily rewritten as 


1 ( d Z H , dh\ 

H VdrF + CJthn d^/’ 


m 


+ yr 


sinh 2 n 0 \ dB 2 


1 d 2 0 


+ cot9 


d©\ m 

d9 / sir 


sin 2 9 


= 0 


( 8 ) 


Since the first and second terms are dependent only on r) and the third and 
fourth terms only on 6 we may reasonably assume that they are independent of 
each other. Therefore each is equal to a constant which by convention is 
chosen to be £(£.+!) 


l/dH 2 . dH\ m 2 

HVdn 2 hn dn) sinh 2 n A(4+1) 

1 /d 2 0 , ft d0\ m 2 _ 

0 [d 9 2 * cotB de)“ sin 2 9 ~ " £( * +1) 

These two equations may be rewritten as 

0 + »«»> f ° 

0 + „ t9 -^]e= o 


f 1 


(9) 


( 10 ) 


(ID 


( 12 ) 
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Equations (11) and (12) may be simplified by the introduction of two new 
variables 

£ = coshn (13 


and 


p e COS0 


(14 


I * 

Considering first equation (11) we see that 


!£ _ 

dn 


sinhn 



coshn 


K 


sinh 2 n = cosh 2 n-l = £ 2 -l 


dH dH d£ , dH 
df| d5 drf sinhn « 


d 2 H 


dry 


d , , . - dH x , d dH 

a? (slnhn) df + s 1 " 1 ’ 1 ’ it,' 3e 


. dH . , d£ d dH 

coshn if + sinhn dTT dTT 


r dH , . , z d 2 H 

^ + Slnh n d £ z 


d 2 H 


dn' 


- (S ! -» M + 5 3? 


d£ : 


dH 

dC 


Substituting equations (15) and (16) into (11) results in 


(? 2 -l) d H 


— v + 2E ^ 

dF + ^ 


- [*(£+1) + ^ ]h 


= 0 


This i« Legendre’s equation with which we will deal later. 


u; 


Cli 


(i 
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Applying the same method to equation (12) we see that 


de " “ sin9 


d 2 p 

c3 2 


cos 0 = - y 


sin 2 0 = 1 - cos z 9 ° 1 - y' 


d 0 d 0 dy . o d 0 

de; = dy d 0 = - sin0 di I 


(18) 


d 2 e . dO . 0 d dO 

— — - e= - coso -rrr - sino -vr 

d0 2 dy d8 dy 


d 0 dy d d 0 

- coso -r— + sino • -r^r i — 3 — 
dy d 6 dy dy 


d0 . 2fl d 2 0 
cos 0 + sin 0 


uV) i!a . u J0 

^ p ; dy 2 p dy 


(19) 


Substituting equations (18) and (19) into equation (12) results in 


0 - f + - 10] 0 ' 0 


( 20 ) 


Multiplying (20) by -1 results in 


o j2 - i > 0 + ^ 


d 0 

dy 


ff - [*« +1 > + ^0] S * 0 


( 21 ) 


Thus the separated equations for H(ri) and 0(0) are of identical form and are 
both Legendre’s equations for which the solutions are associated Legendre 
functions tff the first kind P^(£) and the second kind Q™(£), 
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Thus the separated equations for U(rj) and 9(0) are of idential form and are 
both Legendre's equations for which the solutions are associated Legendre 
functions of the first kind P? (£) and the second kind Q^f)- To solve 
equation (21) we sec m=o and solve for the Eero order Legendre function. 
The first few of which are listed here: 


Po(0 = 1 


Pi(f> = f 


Pa(0 = t,(3C*-l) 


Q 0 (e) = I log [fri] 


Qi(o = eco(e) - i 


q.2 (0 = p 2 («Qo(o - c: e 
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The Legendre functions of order m greater than zero are the Legendre asso- 
ciated functions which are derived from the zero order function by the 
following relationships 


p"(0 


(£ 2 - i >“ /3 


(28) 


and 


Q« (O 


{^ Z _ l ) n /2 


d B Qe(Q 

d£° 


(29) 


These functions are very useful for calculating the electric and magnetic 
fields in and around many objects which are most conveniently described in 
spheroidal coordinates. 
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2.3.6 Electric Scalar Potential Around A Solid Permeable Prolate Spheroid 
Immersed In A Uniform Field - 


2. 3. 6.1 Introduction - 

2. 3. 6. 2 Matching Boundary Conditions At The Surface Of The Spheroid ' 

It is shown in Section 2,3.5 that when Laplace's equation is solved in pro- 
late spheroidal coordinates by the method of separation of variables the 
resulting ordinary differential equations in the variables £ = coshn and 
y = cos9 are Legendre equations. It is further shown there that the eigen- 
functions of Legendre's equation are the Legendre functions of the first 
kind P^COj the Legendre functions of the second kind Q^(£) and the associated 
Legendre functions of the first and second kind P™(£) and Q™(£). The poten- 
tial around a prolate spheroid polarized parallel to its major axis is 
axially symmetric. Therefore only the Legendre functions of the first and 
second kinds are necessary to construct a solution. 


2, 3. 6. 2 Matching Boundary Conditions at the Surface of the Spheroid 

Two potential 'functions, <f>^ inside the spheroid and <f >2 exterior to the 
spheroid, must be constructed. All of the Legendre functions of the second 
kind include the terms 

V ?> = \ log (1) 

or 

V“> “ i ^ (i^) «> 


Mote that both functions in (1) and (2) approach 00 as their arguments 
approach unity. Since the z-axis through the spheroid corresponds to y=l, 
the Q^(u) functions must be excluded from both and Similarly, the 

z-axis inside the spheroid corresponds to £=1, where 


2c 


(3) 
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and therefore Q 2 (£) must be excluded from 

Now that the pathological functions have been excluded we are ready to 
construct well behaved potential functions for cf> ^ and <f > 2 • Using all of the 
allowed eigenfunctions of the axially symmetric Legendre equation we see 
that the general solution for and <}> 2 is 

00 

= z A 5 P„(C)Po(y) (4) 

1 i= o * % % 


” [. 
=0 L 






P£(U) 


(5) 


where 


5 


coshn = 


r l+r 2 

2c 


( 6 ) 


r — r 

y = cos0 = - (7) 

and r^ and r 2 are the magnitudes of the vectors r^ and r 2 from the two foci of 
the spheroid to any point in space, c is the semi focal length of the spheroid, 
and n and 0 are the coordinates used to define the ellipsoidal coordinate sys- 
tem as shown in Section 2.3.2, 

The unknown coefficients A^, and will be evaluated by imposing three 
boundary conditions: 

1) At great distance from the spheroid the potential <J> 2 must approach the 
potential of a uniform field. Since a potential function may have a constant 
added without changing the associated field E because 


£ «= - Vij> 


( 8 ) 


we are free to choose the origin of the spheroidal coordinate system as the 
point of zero potential. The potential function which gives rise to a uniform 
field E 0 in the z direction is 

C9) 


$' = -zE 0 
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But from the transformation from prolate spheroidal to Cartesian coordinates 
derived in Section 1, we see that equation (32) of Section 2.3.2 is 


2 = 

(10) 

- cP 1 (€)P 1 (U) 

(11) 

At great distances from the origin, Q^(£) becomes negligibly 
(5) simplifies to 

small therefore 

l 

Ol> 

(12) 


It is therefore clear that 



B l = -cE 0 (13) 

= 0 for all £>1 


2) At the interface between regions 1 and 2 (the surface of the spheroid 
where £ = £q) the potentials must be equal. From (4), (5) and (13) we see 
that 


OO 00 

WW’ “ - cE „ p i<> i >5o 


(14) 


Since the Legendre Polynomials P^(y) ate orthogonal over the interval -1 <_ y <_ 1 
we can multiply (14) by a Legendre polynomial P^OO and integrate from -1 to +1 
to separate the coefficients. This operation results in 

« 

A oW - c oW ■ 0 < 15 > 


«•% 
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A 1^0 ■ C 1^1^0* = _cE °^o 

(16) 

A k p k (C o ) ~ c k W 0 0 for a11 k>1 

(17) 


3) At the interface between the two regions 1 and 2 the normal derivatives 

of the potentials multiplied by the relative permeativitics of their respective 

regions must be e'jval. This is required to insure continuity of the normal 
component of flux across the boundary. Thus both Dirichlet and Neumann 
boundary conditions are imposed. Thus 

94> 94) 

e r = IT Wh£n ? “ K 0 < 1S > 

or from (A) and (5) 




E 

s,=0 


D 


B £, P f/V 


+ C 


*«o>] 


p £ (y) 


(19) 


Again invoking the orthogonality of the Legendre polynomials over the interval 
~l£p£l and using (13) we can separate the coefficients with the result: | 


Er A 0 P 0''V “ C 0%^(P ~ 0 ( 2 °) 

e r A 1 - B^^Cq) - C lQ £(C 0 ) = 0 (21) 

t 

e r A k P k^0' ) “ C k Q k^0 } = 0 (22) 

Equations (20) and (21) can be simplified by recalling that 

P 0 (O - 1 P'(D = 0 (23) 

P x (?) = € Pj(€) = 1 


(24) 
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Since Pq(O“ 0, we conclude that 


From (15) we see that Aq“0 


(2i 


V 



At this point we shall deal with all of the coefficients for k>l recalling 
(17) and (22) 


W 5 0> - °kW " 0 

Vk <5 0> - c k Q k <5 0> ■ 0 

Multiplying (17) by "Q^^q) and (22) by Q^C^g) and adding the two equaclons, 
we find that 

\ - p k< e o^.;«oD - ° C2. 


(17) 

( 22 ) 


If we can prove that the term in the brackets is never equal to zero we have 

The Russian 


established that all and coefficients are zero for k>l. 


mathematician N.N. Lebedev 1 shows that the Wronskian W(P^(£) ,Q^(£) ) which is 


defined as 


W 





P £ (?) Q £ «> 

pjce> qj(o 

p a (C)Q £(C) - p£(OQ £ (« 


is equal to l/(l-£ 2 ). Therefore, equation (26) can be rewritten as 




(2 


(2 
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Now £.>1 and the Q. functions are always positive for £>1. Furthermore, all 

U K 

of the roots of He in the interval -1 <. £ <_ 1. Thus P k (£ Q ) >0 

Cq^I. Since the term enclosed in the brackets Ib always greater than zero 
the inly possible solution to (28) is 




=■ 0 for ell k>l 


(29) 


from which we deduce from (17) that 


C k ° 0 for all k>l. 


(30) 


Thus we see that the only non-zero coefficients of (4) and (5) are A^, 
and C^. But it has already been shown by the first boundary condition that 
B.=-cEn Therefore, equation (21) simplifies to 


Gy A^ — C^Q^(5 q) 1=1 ~nEo 


(31) 


But 


d£ Q x (0 = Q 0 U) - 


(32) 


Substituting (32) into (31) results in 
E r A i " c i [q 0 ~ _ cE o 




(33) 


or 


°i r £ "I cE ° 

A i ■ * [«o - ffci] — 




(34) 


Substituting (34) into (16) 


r T [jio ' t^i] - C A 


cE 0 C 

— < e *" l > * ^0 


(35) 


S3SSB5T«T^ 






But 
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Ql'(C) - %(Q -1 


1‘herefore (35) simplifies to 


C 1 |j e r rl )Q 1 + “ ^cEo(e r -l) 




(36) 


Dividing (36) by (e r -l) results in 


c?E, 


9 i (0 + fe 


. e«e 


0 


(37) 


Substituting (37) into (16) results in 


c ^oQi 


A x = - cE 0 + 


Q 1 + ( er -im z -l) 


* . 


(38) 


which simplifies to 


- cE c 


1 1 + (e r -l)(5 2 -l)Q l (€) 




(39) 


Thus we have evaluated and in terms of c, E 0 , e r and £q» and we can write 
the solution for the potential in the two regions 1 and 2 as 


$1 = 




$ 2 = c i.Qi’<Ou - cEoCh , ? 0 <C<” 


(40) 

(41) 


But by the transformation which generated the prolate spheroidal coordinates 
resulted in the relationship 


z = c£p (1-35) 

Therefore we see that the field within the spheroid which is given by 

t s • V<£ 


is uniform and is given by 

A, 

z 


p 1 • 

C SI 7 . 


(42) 
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10)1 ‘ 


Wo have solved the case of a dielectric spheroid immersed in an otherwise 
uniform electric field. To specialize to the case of a conducting wire we 
allow the dielectric constant to go to infinity. Under this condition Ai=0 


c£qEq 

Qi(S*) 


We substitute (13) and (43) into (4) with the result 


= =Eo 


l site.) { J " 


Recalling (12) , (13) , (18) and (19) from Section 2.3.2 


n = yx 2 + y 2 + (z-c) 


rj = yx 2 + y 2 + (z+c) 


ri + r 2 
2c 


- r i 
2c 







Page 109 


we see that for every point x, y, z in space We can calculate the electric 
potential <f>. This equation is used to determine the boundary conditions 
for the numerical solution of the electric field around the wire in the 
presence of plasma. A plot of the electric field around a conducting 
spheroid is shown in Figure 2.3.3, The electric flux lines were calculated 
from an analog of the magnetic vector potential 


D = VxA, (45) 

A smaller plot for a thin wire is shown in Figure 2.3.4. Equation (44) is 
also used to initialize all of the interior grid points before the itera- 
tive solution is begun. The details of this variable grid iterative solu- 
tion will be described in the next report. 


t 
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Figure 2.3.3 Electric field around a conducting 
spheroid. 
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2.4 Potential Build-Up Due To Plasma Contactor Failure 

In an easterly orbit, with the satellite deployed upward, the satel- 
lite collects electrons and the Shuttle collects ions. Because of the low 
thermal velocity of the ions, it is necessary to have ■* plasma contactor on 
board the Shuttle in order to maintain it at plasma potential. If the 

9 

plasma contactor fails, the Shuttle will acquire a negative potential as a 
result of the inability of the Orbiter to neutralize the electrons coming 
down the wire from the satellite. If the satellite has a plasma contactor 
it will be a better collector and the Shuttle will acquire a large poten- 
tial . 

Two simulations have been run in order to estimate the potential ac- 
quired by the Shuttle if its plasma contactor fails. These simulations are 
relevant to the first electrodynamic mission with the Shuttle on an orbit 
| of 28.5° inclination and 295 km altitude. The satellite is deployed upward 

on a 20 km long tether. These simulations refine and confirm the value of 
-2kvolt for the Shuttle potential as it was estimated in the Quarterly 
Report #1. In the first one the satellite has no plasma contactor and is 
modeled as a 1.5 meter disameter sphere. The Shuttle is modeled as a 12.4 
meter diameter sphere. The wire is 20 km long and has a resistance of 4000 
ohms. The Shuttle acquired a voltage of about -1195 volts and the satel- 
lite was at +1206 volts. The wire current was .253 amps. In the second 
simulation, the satellite has a plasma contactor and the other parameters 
| are the same. The Shuttle acquired a potential of -2064 volts and the 

| satellite was within a volt of plasma potential. The current in the wire 
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3.0 PROBLEMS ENCOUNTERED DURINC REPORTING PERIOD 


None . 


4.0 ACTIVITY PLANNED FOR THE NEXT REPORTING PERIOD 


During the next reporting period we will carry out the software modi- 
fication of the force and torque package in our rotational dynamics com- 
puter code. The implementation of the reeling control laws will then allow 
us to simulate actual maneuvers with a good fidelity. Our plotting program 
will also be modified because presently it is not formatted to receive and 
Interpret the rotational dynamics variables. 

The three-dimensional slack tether computer code SLACK3 is now avail- 
able and tested. During the next reporting period it will be used to 
simulate actual cases of tether break and to investigate suitable Shuttle 
avoidance maneuvers. In the mean time the analytic studies on the slack 
tether will continue. The major goal is the evaluation of the wrinkle’s 
typical scale size at the onset of the compressional instability. 

The development of the computer code to evaluate the electric field 
around the severed tether is in progress. The next reporting period will 
be devoted to the implementation of the electric potential variations due 
to the plasma. Presently the program can evaluate the electric potential 
around a charged body in vacuo . 


